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We  discuss  an  algorithm  for  computing  the  stationary 
probability  vector  of  an  infinite-state  Markov  chain  whose 
transition  probability  matrix  has  a block-partitioned  struc- 
ture. Such  matrices  arise  in  a wide  variety  of  queueing 
models  as  well  as  generalized  random  walk  problems.  Tradi- 
tionally, the  analytic  approach  to  this  type  of  problem  has 
been  through  complex  variable  methods.  We  present  an  alter- 
nate and  unified  treatment  of  this  problem  and  obtain  an 
algorithm  which  utilizes  only  real  arithmetic  computations. 

In  addition,  most  of  the  intermediate  steps  of  the  algorithm 
have  useful  probabilistic  interpretations. 

We  obtain  an  adequate  number  of  the  initial  compon- 
ents of  the  invariant  vector  by  using  a purely  probabilistic 
argument.  Higher  components  are  evaluated  by  matrix-itera- 
tive methods.  The  first  and  second  moments  of  the  stationary 
distribution  are  also  found  in  computationally  tractable 
forms.  The  APL  program  used  to  implement  the  algorithm  is 
listed  and  several  numerical  examples  are  presented. 
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GLOSSARY 


A - the  irreducible  stochastic  matrix  V A,, 

u „ v 
v = 0 


{An)n=0 


A*  (z) 


A*(z) 


A 1 
Ak 


{Bv} v=0 


V 


Ej 


a sequence  of  nonnegative  mxm  matrices  whose  sum 
is  stochastic  and  irreducible. 


the  matrix  generating  function  J A zv. 

v=0 

the  matrix  obtained  by  differentiating  the  matrix 
A*(z)  entrywise,  v times  with  respect  to  z. 


- Av  (I-A-,  ) 


-1 


a sequence  of  nonnegative  matrices  satisfying 

co 

B0®  + I Bye  = e where  Bq  is  n*n  and  is  nxm 
v=l 

for  k>l. 


(*0Bk  + *iV  (I-A1> 


-1 


the  invariant  probability  vector  of  the  matrix  L. 
the  mean  vector 


- 

dz  z=l~  . 

J 


- a vector  (of  appropriate  dimension)  each  of  whose 
components  is  1. 

- the  mean  recurrence  time  for  the  state  (0,j). 


GOO 


G(z) 


G 

a 

H 

H (k) 


H(z) 


h* 


I 


- the  matrix  { G j j , (k)  } where  Gjjt  (k)  equals  the 
probability  that  starting  in  state  (i+l,j),  ilO, 
the  level  i is  reached  for  the  first  time  ir 
state  ( i , j * ) in  exactly  k transitions. 

00 

. r k . 

- the  transform  matrix  \ G(k)z  satisfying  the 

k=0  °° 

functional  equation  G(z)  = £ zAvGv(z). 

v=0 

- the  matrix  {Gjji}  where  Gjji=9j*»  for  11  j,  j'lm. 

- the  stationary  probability  vector  of  G. 


- H ( z) 


z=l* 


- the  matrix  { H ^ ^ i (k) } , where  Hjj,(k)  equals  the 

probability  that  starting  in  state  (l,j),  the 
level  0^  is  reached  for  the  first  time  in  state 
(0 , j 1 ) in  exactly  k transitions. 


n 

- the  transform  matrix  i H(k)z  . 

k=0 


- the  mean  vector 


- the  mxm  identity  matrix. 


i - the  "level  i",  where  ill,  consisting  of  the  set  of 

states  {(i,j),  lljlm}  in  the  infinite  Markov 
chain  P. 


■ 


K (k) 


K ( z) 


L(k) 


L(z) 


X(z) 


X(n)  (z) 


u(z) 


u(n) (z) 


v(z) 


the  matrix  {Kjj,(k)}  where  Kj^,(k)  equals  the 

probability  that  starting  in  state  (l,j),  the 
Markov  chain  returns  to  level  1 for  the  first 
time  in  the  state  ( 1 , j * ) in  exactly  k transi- 
tions . 


the  transform  matrix  l K(k)z  . 

k=0 


- n*>i2=r 


{L j j , (k)  } where  Ljj,(k)  equals  the  probability 

that  starting  in  state  (0,j),  the  Markov  chain 
returns  to  level  ()  for  the  first  time  in  the 
state  ( 0 , j ' ) in  exactly  k transitions. 


the  transform  matrix  £ L(k)z^. 

k=0 

the  transition  probability  matrix  of  the  infinite 
Markov  chain  having  the  particular  structure  of 
interest. 


co 

the  vector  generating  function  £ xvzv. 

v=l 

the  vector  obtained  by  differentiating  the  vector 
X(z),  n times  with  respect  to  z. 


- the  invariant  probability  vector  of  the  matrix  P. 


the  appropriately  normalized  right  eigenvector 
of  the  matrix  A*(z),  corresponding  to  the  Per.ron- 
Frobenius  eigenvalue. 


the  vector  obtained  by  differentiating  the  vector 
u(z) , n times  with  respect  to  z. 


the  appropriately  normalized  left  eigenvector  of 
the  matrix  A*(z),  corresponding  to  the  Perron- 
Frobenius  eigenvalue. 


i 


v<n)  (z) 


the  vector  obtained  by  differentiating  the 
vector  v(z)  , n times  with  respect  to  z. 


1 vAv®.  ■ 
v=0 


- the  Perron-Frobenius  eigenvalue  of  the  matrix 
A* (z)  . 


' (z)  - the  n-th  derivative  of  6(z)  with  respect  to  z, 


- a diagonal  matrix  with  the  elements  of  x along 
the  diagonal. 


- the  invariant  probability  vector  of  the  matrix  K. 


- the  mean  vector  K(z)  , e. 

dz  z=l— 


- the  mean  vector  ^ G(z)  z_^e 


- the  invariant  probability  vector  of  the  matrix  A. 


- the  matrix  { II  . . , } where  IU.it  = tt  • f 


- the  quantity  ti_6, 


1.  INTRODUCTION 

We  are  concerned  with  a class  of  infinite  Markov 
chains  with  stationary  transition  probabilities,  having  a 


transition  matrix  P of  the  following  form: 


where  the  matrices  Av , v>0,  are  square  substochastic 

matrices  of  order  m.  The  matrices  Bv,  v>l,  are  nxm,  while 

Bq  is  nxn  and  Cg  is  mxn.  The  state  space  of  this  Markov 

chain  is  the  set  {(0,j),  lijin  and  (i,j),  i£l,  l£j£m}. 

00 

The  matrix  A = £ Av  is  stochastic  and 

, v=0 

j. 

CO 

(2)  C0e  + l Aye  = e, 

V = 1 

oo 

Bo®  + l 3v£  = 

V = 1 

where  e = (1,1,...,!)'.  We  shall  derive  an  algorithm  to 


1 


f 


i 


2 


compute  the  steady-state  probability  vector  x = 

(xq , xi ,x2 / • . • ) of  the  matrix  P,  where  Xg  is  an  n-vector 
and  x^,  k>l,  are  m-vectors.  This  amounts  to  solving  the 
infinite  system  of  linear  equations 

(3)  x P = x,  x e = 1 , 

or  equivalently, 

<4>  *0  = *0B0  + *1C0 

k+1 

*k  = x0Bk  + l xvAk+1_v,  for  k>l. 
v = l 

We  define  the  probability  generating  function  vector  X(z) , 
O^z^l,  as 

00  k 

(5)  X ( z)  = l xkz 

k=l 

» 00  k+1 

= *0  l Bkzk  + l zk  l *vAk+l-v, 

k=l  k=l  v=l 

and  the  generating  function  A*(z),  Oizil,  of  the  sequence 
of  matrices  { An } as 

Qb 

(6)  A* (z)  = l Avzv . 

v = 0 

If  we  interchange  the  order  of  summation  in  the  second 
term  on  the  right  of  Equation  (5) , we  have 


It  is  traditional  to  attempt  to  derive  the  vector  Xq  from 
Equation  (7)  by  using  complex  variable  methods,  based  on 
an  application  of  Rouche's  theorem.  In  practice,  however, 
this  method  may  lead  to  highly  unstable  numerical  computa- 
tions. We  shall  derive  the  vector  Xq  using  a purely  pro- 
babilistic argument.  It  should  be  stressed  that  our 
approach  will  utilize  only  real-arithmetic  algorithms  and 
so  avoids  many  of  the  numerical  problems  associated  with 
the  complex  variable  methodology.  Our  discussion  reviews 
and  generalizes  a number  of  earlier  results,  used  in  the 
analysis  of  specific  queueing  models  [24,25,27]. 

Markov  chains  of  the  type  (1)  appear  as  the  em- 
bedded Markov  chains  in  a large  number  of  queueing  models. 
Computing  the  vector  x is  a crucial  step  in  the  numerical 
evaluation  of  many  quantities  and  probability  distribu- 
tions of  relevance  to  the  theory  of  queues.  A list  of  sub- 
stantially different  queueing  models,  which  are  amenable  to 
the  present  analysis  is  given  in  Section  8. 


I 


2.  THE  FIRST  PASSAGE  TIMES  FROM  LEVEL  i+1  TO  LEVEL  i 

Consider  the  first  passage  times  from  the  set  of 
states  i+1  = { (i+1, j) , je (1, . . . ,m) } to  the  set  of  states 
i.  = {(i#j)»je(l#...»m)},  i-1 . The  set  i^  will  henceforth 
be  referred  to  as  level  _i.  Let  G j j i 00  , 1- j » j 1 -m,  k>l,  be 
the  conditional  probability  that,  starting  in  the  state 
(i+l,j),  the  process  reaches  the  level  i for  the  first 
time  in  the  state  ( i , j ' ) after  exactly  k transitions. 

Define  the  sequence  of  matrices  {G(k),k£l}  such  that  G(k)= 
(Gjji(k)}.  This  sequence  of  matrices  defines  completely 
the  first  passage  time  distributions  from  states  in  the 
level  i+1  to  the  level  _i.  These  matrices  were  studied  in 
great  detail  by  Neuts [21 , 24 ] . We  now  review  a number  of  im- 
portant results  from  these  papers,  which  are  needed  in  the 
sequel . 

The  matrices  G(k),k_l,  are  most  conveniently  stu- 
died by  considering  the  matrix  of  transforms  G(z),  defined 
by 

CD 

(8)  G (z)  = l G(k)zk,  for  Oiz<l. 

k=l 

By  using  a standard  first  passage  argument,  it  is  shown 
that  G ( z)  satisfies  the  matrix  functional  equation 


4 


5 


(9)  G ( z)  = l zA  G V ( z ) . 

v=0 

In  [21]  it  is  shown  that  Equation  (9)  uniquely 
determines  the  sequence  of  matrices  (G(k) ,k-l}.  For  the 
process  eventually  to  reach  level  _i  from  any  state  in  i+1 , 
the  matrix  G=G(1-)  must  be  stochastic.  The  following 
theorem  is  proved  in  [21] , assuming  certain  irreducibility 
conditions,  which  are  generally  satisfied  in  applications 
and  which  we  shall  not  repeat  here. 

Theorem  1;  Let  n_  be  the  invariant  probability  vector  of 
the  irreducible  stochastic  matrix  A,  i.e.,  the  unique  solu- 
tion to  the  equations 

(10)  £ A = JL  and  £ e = 1. 

oo 

Also  let  £ = [ vAve.  Then  the  equation 

v=l 

oo 

(11)  G = l AVGV, 

v=0 

has  a minimal  nonnegative  solution  which  is  stochastic  if 
and  only  if  p = £ §_  -1.  The  matrix  G is  then  also  the 
unique  nonnegative  matrix  satisfying  that  equation. 

Remarks : 

a)  For  what  follows,  we  assume  the  above  mentic  ed 
irreducibility  conditions  hold. 
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b)  In  the  queueing  context,  p as  defined  above, 
represents  the  expected  number  of  arrivals  during  a suita- 
bly averaged  service  time  and  corresponds  to  the  familiar 
traffic  intensity . 

c)  If  p*l,  the  chain  is  null-recurrent  or  transi- 
ent and  therefore  no  solution  to  the  equations  (3)  exists. 
In  the  sequel,  only  the  case  p<l  is  discussed. 

d)  The  matrix  G may  most  conveniently  be  computed 
by  modified  successive  substitutions.  This  corresponds  to 
successively  evaluating  the  matrices 

(12)  G(0)  = (I-A1)_1A0, 

00  , 

G (k+1)  = l ( I-A,  ) -avg'j  (k)  , for  k>0. 
v=0 

v^l 

It  was  shown  in  [21],  that  this  sequence  is  entry-wise 
strictly  increasing  and  converges  to  the  matrix  G. 

For  future  reference,  we  introduce  the  vector  2 of 
stationary  probabilities  corresponding  to  the  stochastic 
matrix  G and  the  square  matrix  G of  order  m,  whose  rows 
are  all  identical  and  equal  to  the  vector  2*  Since  the 
irreducible  matrix  G has  spectral  radius  equal  to  one,  the 
matrix  (I-G)  is  singular.  It  is  an  important  and  well- 
known  result  though,  that  the  matrix  (I-G+G)  is  non-singu- 
lar. (see  Kemeny  and  Snell  [11]).  We  shall  also  need  the 
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mean  vector  £ defined  by 


(13)  £ = l k G (k)  e . 

k=l 


Theorem  2 : The  vector  £ is  given  by 
(14)  v = (I-G+G) [I-A+G-A (B)G]_1e 

where  A (£)  is  a diagonal  matrix  of  order  m,  with  diagonal 
entries  elf  e2,  •••,  Bm. 


Proof : 


(15) 


~ [dJ  G(z)J  z-iS 

~ \ l AvG  (z^+  I zAv  l Gr  (z)  G'  (z)GV-r_1(z)l__-,e 
L v=0  v=l  r=0  Jz-i 

= \g  + l A Y GrG'  (1)  GV-r-1~|e 
L v=l  r=0 


<*>  v-1 

£ + I Av  I Grw. 
v=l  r=0 


<*>  v-1  oo 

Now  I Au  I Gr  ( I-G+G)  = l A fI-GV+vG) 

1 n i v 


v=l  r=0 


v=l 


= A-G+A (e)G. 


Therefore 


(16)  (I  - (A-G+A (£)G) (I-G+G) _1] u = e, 
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and  after  writing  I = (I-G+G)  ( I-G+G ) "'’and  simplifying,  the 
desired  result  follows.  Neuts  has  shown  in  [24],  Theorem 
4,  that  the  inverse  used  in  Formula  (14)  exists.  Note  that 
that  Pj  equals  the  expected  number  of  transitions  during  a 
first  passage  from  the  state  (i+l,j)  to  the  level  i. 


Corollary  1:  The  inner  product  £ £ is  given  by 
U7)  £ £ - x^-  • 

Proof : Since  (l-p)g  y = ( 1-p ) £[ I-A+G-A (£) G] -1e  and 

£[  I-A+G-A  (8)  G]  = (l-p)g,  it  follows  that  (l-p)£  £ = £ e = 1* 

Corollary  1 provides  a powerful  accuracy  check  on 
our  numerical  computations  as  well  as  having  its  own  pro- 
babilistic significance,  i.e.,  £ £ equals  the  "average" 
number  of  transitions  required  to  go  from  level  i-H  to 


level  i. 


3.  THE  FIRST  PASSAGE  DISTRIBUTIONS  FOR  LEVELS  0 AND  1 


Let  us  define  Ljj' (k)  to  be  the  conditional  proba- 
bility that  starting  in  state  (0,j),  the  process  returns 
to  level  0 for  the  first  time  in  state  ( 0 , j ' ) after 
exactly  k transitions.  Let  {L(k),k;Ll}  be  the  sequence  of 
matrices  L (k) ={L j j i (k) } , which  completely  defines  the  dis- 
tribution of  the  first  passage  times  from  level  ()  back  to 
level  0.  We  analogously  define  the  sequences  of  matrices 
{H(k),k£l}  and  {K(k),k>l}  as  the  first  passage  distribu- 
tions from  level  1 to  level  0^  and  from  level  1^  to  level  1, 
respectively.  In  the  context  of  queueing  theory  {H(k),k£l} 
corresponds  to  the  densities  of  the  number  of  customers 
served  during  busy  periods  starting  with  one  customer.  We 
define  the  corresponding  matrix  generating  functions  as 
follows : 

CO  00 

(18)  L (z)  = l L (k) zk , H ( z)  = l H(k)zk, 
k=l  k=l 

oo 

K ( z)  = l K(k)zk,  for  0<z<l. 
k=l 

It  follows  by  standard  first  entrance  arguments, 
that  the  following  equations  hold: 


9 


10 


(19)  L(z)  = 


zBg  + l zB^G V 1 (z)H(z) , 
v = l 


i 


-1 

H ( z)  = zCq  + l zA^G V (z)H(z) 
v=l 


1 

z[I  - I zAvGv_1(z)]  CQ/ 
v = l 


.v-1 


K ( z)  = zCq  l zrBj  l zBvG  (z)  + l zA  Gv_i(z) 
r=0  v=l  v=l 


-1  °°  -1  °°  -1 
= zCg (I-zBg)  l zBvGV~  (z)  + l zAvGV"  (z) 

V=1  V=1 


To  show  that  the  inverse  in  (19)  exists,  we  see  that 

OO  00 

£ zAvGV  ■'"(z)  e K l A GV_'*'e  < (A-Ag)e.  But  under  the 
v=l  v=l 

assumed  irreducibility  conditions,  which  normally  hold  in 
practice,  the  matrix  (A-Aq)  is  strictly  substochastic.  By 

Corollary  2.2  in  the  appendix  in  Karlin  and  Taylor  [10] , we 

00  - 

have  that  the  matrix  £ zA  Gv-i(z),  0<z<l,  has  spectral 

v=l 

radius  less  than  one  and  therefore  the  desired  inverse 
exists.  Note  that  the  matrices  L=L(1-),  H=H(1-),  K=K(1-) 
are  all  stochastic,  whenever  G is  stochastic.  For  example 
this  is  verified  for  H as  follows.  Since  CQe  = Age,  the 
vector  He  may  be  written  as 

°o  - 

He  = (I  - l AVGV'1) " A0e, 

V = 1 
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oo 


but  clearly  Equation  (11)  implies  that  G=(I- 
and  therefore  He  = Ge  = e. 


v 


-1,-1. 

) Ar 


We  define  the  invariant  probability  vectors  d and 
of  the  matrices  L and  K as  follows: 

(20)  dL(l)  = d and  d e = 1 , 

<K(1)  = < and  < e = 1. 


In  the  sequel,  we  shall  derive  explicit  expressions  for 
the  mean  vectors  d* , h*  and  k*  defined  by 

(21)  d*  = L'(l)e,  h*  = H ' ( 1 ) e , k*  = K'(l)e. 


Theorem  3:  Provided  the  vector  £ vB  e is  finite,  the  mean 

v=l 

vectors  h* , d*  and  are  given  by 


(22)  h*  = 


= U-  l AvGv  1)"1(e+f  l A - l A Gv-1+  l (v-l)Ava| 
v=l  [ [y=l  V v=l  v=2  J 


(I-G+G) "1u 


> » 


00  - oo  OO  to  ^ 

= £+  l BvGV-1h*+  l Bv-  l BVGV_1+  l (V-1)B  G 
V = 1 1_V=1  v = l V = 2 


(I-G+G) 


-1 


P , 


and 


K*  = e + CQ (I-Bq) + C0(I-B0)-1f  l Bv-  l BVGV-1 

l Lv=i  v=1 

+ l (v-l)BvG I + l Av-  l AVGV_1 
v-2  J v = l v=l 


+ J ( v-1) AG  > (I-G+G)~V 

v=2 


Proof : If  we  differentiate  H(z)  with  respect  to  z we 


obtain 


_ oo 

H'(z)  = CQ  + l AvGV_1(z)H(z)  + l zAvGV"1(z)H' (z) 


» v-2 

+ I zAv  I Gr(z)G' (z)GV"r“2(z)H(z) . 
v=2  r=0 


Letting  z tend  to  1-  yields 

H ' ( 1)  = H + f a,GV'1  H'(l)  + I Av  l GrG'  (l)GV_r"2H. 
V=1  v=2  r=0 


Therefore  we  have 


(23)  h*  = H 1 ( 1) e = e + l AvGV_1h*  + £ Av  £ Grp. 

v=l  v=2  r=0 


00  V “ 2 00 

1 Av  l Gr(I-G+G)  = l A (I-GV-1+ ( v-1) G) 

v=2  r=0  v=2 


00  oo 


“ l l AVG  “ + l (v-l)AvG. 

v=l  v=l  v=2 
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We  see  that  the  formula  for  jc*  involves  the  inverse 
of  the  matrix  (I-Bg) . If  the  matrix  (I-Bg)  were  singular 
there  would  exist  a relabeling  of  the  rows  and  columns  of 
Bq , such  that  it  may  be  written  in  the  form 


Clearly,  a subset  of  the  states  (0 , 1) , . . . , (0 ,n)  would  form 
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an  irreducible  class  and  the  infinite  matrix  P would  then 
be  reducible.  In  the  irreducible  case  under  consideration, 
the  matrix  (I-Bq)  is  necessarily  nonsingular. 


1 


. THE  STATIONARY  PROBABILITY  VECTOR  OF  THE  MATRIX  P 
Let  x(i,j)  be  the  limiting  probability  that 
immediately  following  a transition,  the  Markov  chain  is  in 
the  state  (i,j).  In  the  positive  recurrent  case  under  dis- 
cussion the  quantities  x(i,j)  form  an  infinite  probability 
vector  x,  which  we  will  write  in  the  partitioned  form  (xq, 

— 1'— 2'-**)*  The  system  of  equations  (3)  has  a unique  solu- 
tion with  all  x(i,j)>0  if  and  only  if  pel. 

Using  classical  arguments  in  the  theory  of  Markov 
renewal  processes,  we  shall  derive  explicit  expressions  for 
the  vectors  Xq  and  x^ . 

Theorem  4 : The  vectors  Xq  and  x-^  are  given  by 
(24)  -0  = dd*  ' -1  = k|*  * 

Proof : If  we  consider  each  transition  in  the  infinite 
Markov  chain  P as  a discrete  time  step,  then  the  times 
between  successive  visits  to  the  states  (0,1) , . . . , (0,n) 
and  the  states  visited,  define  a Markov  renewal  process 
with  n states.  The  sojourn  times  in  this  n-state  Markov 
renewal  process  are  lattice  random  variables  and  the  tran- 
sition probability  matrix  of  the  process  is  given  in  an 
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equivalent  form  by  the  matrix  of  generating  functions  L(z). 
There  is  a classical  theorem,  see  e.g.  Hunter  [9],  p.  196, 
in  the  theory  of  Markov  renewal  processes  that  states  that 
the  mean  recurrence  time  Ej  of  a particular  state  (0,j)  is 
given  by 

i n 

(25)  E.  = l d d*  , for  l<j<n, 

i d-;  L,  v v J ' 

J J v=l 

where  dv  and  d*  are  the  v-th  components  of  the  vectors  d 
and  d*  respectively. 

The  mean  recurrence  time  Ej  in  this  finite  state 
Markov  renewal  process  is  none  other  than  the  expected  num- 
ber of  transitions  between  successive  returns  to  the  state 
(0,j)  in  the  infinite  state  Markov  chain  P.  Clearly,  we 
see  that  the  stationary  probability  x(0,j)  is  given  by 


(26) 


x ( 0 , j ) 


for  l<j<n. 


or  in  vector  notation, 
d 

-0  “ dd*  • 

To  derive  an  explicit  expression  for  x-^ , we  consider  the 
Markov  renewal  process  defined  by  the  times  between  suc- 
cessive visits  to  the  states  (1 , 1) , . . . , (1 ,m)  and  the  states 
visited.  The  transition  probability  matrix  for  this  process 
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is  given  in  an  equivalent  form  by  the  matrix  of  generating 
functions,  K(z).  Using  a completely  analogous  argument, 
we  see  that  Equation  (24)  for  x-^  holds. 

Since  the  matrices  L and  K are  known  in  terms  of 
matrices  which  may  be  computed  explicitly.  Formula  (24) 
yields  the  vectors  Xq  and  x-^  in  a tractable  form.  Knowing 
that  the  vector  x must  satisfy  Equation  (3),  we  must  show 
that  the  following  relationship  between  Xq  and  x-^  hold. 

Corollary  2:  The  vector  Xq  is  related  to  the  vector  x^  by 
the  following  equality: 

(27)  xQ  = XqBq  + XjCq 

= XiCod-Bo)-1. 

Proof : By  a lengthy,  but  straightforward  calculation  given 
in  Appendix  I. 

We  see  that  Equation  (27)  provides  us  with  yet 
another  accuracy  check  on  our  numerical  computations. 


5.  THE  DERIVATIVES  OF  THE  PERRON-FROBENIUS  EIGENVALUE 


In  deriving  the  moments  of  the  stationary  distribu- 
tion, we  will  need  explicit  expressions  for  the  derivatives 
of  the  Perron-Frobenius  eigenvalues  and  the  associated 
eigenvectors  of  the  matrix  A*(z),  defined  in  (6).  In  this 
section,  we  derive  the  necessary  recurrence  relations  need- 
ed in  the  computation  of  these  derivatives. 

For  zil,  the  matrix  A*(z)  has  a uniquely  defined 
Perron-Frobenius  eigenvalue  <S(z).  Let  u(z)  and  v(z)  be 
the  corresponding  right  and  left  eigenvectors,  respectively, 
such  that  the  normalizing  conditions 

(28)  v(z)u(z)  = v(z)e  = 1, 

v ( 1)  = it,  and  u ( 1)  = e , 

hold  in  addition  to  the  defining  relations 

(29)  [A*  (z)  -6  (z)  I]  u (z)  - v (z)  [A*  (z) -6  (z)  I]  =0. 

We  denote  by  A*(z),  the  matrix  obtained  by  differentiating 
each  entry  of  A*(z),  v times. 

Theorem  5 : The  derivatives  6 ^ (1)  , u^  (1)  , v^  (1)  , niO, 
may  be  computed  recursively  for  each  n for  which  A*(l)  is 
finite.  The  recursion  formulas  are 
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(30)  6 (0>  (1)  = 1,  u(0)(l)  = e,  v(0,(l)  = tt 


6 d)  = itA^  (1)  e = 


u(1,(l)  = (I-A+n)'1[A*(l)-6(1) (l)I]e 


v(1) (1) 


and  for  n£2 
(31)  6 (n)  (1) 


= (I-A+n)  U - pe. 


= ii_[A*  (l)-6  (1)  (1)  I]  (I-A+JI)  1 
= rrA*  (1)  (I-A+ll)  _1  - pn. 


n (n-v) 

l (v)iAC(1>H.  (D 

v=l 


nrx  (n-v)  (v) 

- I ITU.  (1)  6 (1) 

v = l 


u(nl(l) 


(n-v) 


(I-A+n)  ^(v)[a*(1)-6(V)  (l)l]u  (ij1 

v = l 


■ "l1  (")v(v)  (Du'ui*  e 

_v=l 


vln) (1) 


l <;>v<vl  (i)fas.v(i)-6  (1)  0 (I-A+n) 

vc0 


where  n is  the  square  matrix  of  order  m,  whose  rows  are  all 
identical  and  are  equal  to  the  vector  it . 


Proof ; The  values  corresponding  to  n=0  are  obvious.  By 
differentiating  n times  in  Formula  (29) , we  obtain 
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n (n-v) 

l ($)  [A*  (z)  -6  v (z)  I]  u (z)  = 0. 

v=0 


Premultiplying  (32)  by  v(z) , letting  z tend  to  1 and 
rearranging  terms  leads  to 


n . . n-1  (n-v)  (..) 

(33)  6 n (1)  = l (")  iLA*(l)u(n-v)-(l)-  l iru  (1)  6 (1) 

v=l  v=l 


Letting  z tend  to  1 in  Equation  (32)  and  rearranging  terms 


leads  to 


(34)  (I-A)u(n)(l)  = l (5)  [a*  (l)-6  (v)  (1)  f]u(l(l)')  , 

v — 1 

which  is  a singular  system  of  equations.  Adding  nu^nMl)  = 
[_n  u(n)  (l)]e  to  both  sides  and  noting  that  (I-A+n)-'1-e  = e, 
we  obtain 


/ % , n , x -i  (n-v) 

(35)  u(n)  (1)  = (I-A+n)  1 l (")(_A*(l)-6(v)  (l)lju(l) 

v = l 

+ ( jl  u(»)  (1)  )e. 

In  order  to  determine  u^n^  (1)  in  terms  of  earlier  terms  of 
the  recurrence,  we  differentiate  n times  in  the  normalizing 
condition  v(z)u(z)=l,  and  let  z tend  to  1 to  obtain 


ir  u<n>  (1)  = 


r ,n,  (n)  (n-v) 

- L U)v(1)u  (i) 

v=l 
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Note  that  for  n=l,  jtu^  (1)  = 0,  and 
(37)  v(n)  (l)u(l)  = v(n)  ( 1 ) e = 0. 


Substitution  of  (36)  into  (35)  yields  the  stated  formula 
for  u(n) (1) . The  vectors  v(n) (1)  are  obtained  by  differ- 
entiating the  second  equation  in  (29)  n times  and  setting 
z=l.  We  get 


(38) 


v(n)(l)(I~A)  = l1  (")v(v)  (1)  [a*  (1)  - 6 l?l),)  I 


v=0 


V — 


(n'-v)  \ 


n-v 


but  since  v(n)(l)n  = [v(n)(l)e]^  = 0,  we  have 


v(n)(l)  - "j1 
v=0 


(n-v)  1 

A*_v  (!)  - 6 (1)  I 


(I-A+n) 


-1 


Setting  n-1,  we  obtain  the  stated  explicit  formulas 


6.  THE  MOMENTS  OF  THE  STATIONARY  DISTRIBUTION 


In  the  next  section,  we  shall  develop  a recursive 
algorithm  to  compute  further  components  of  the  vector  x. 

As  a criterion  for  truncation  in  the  infinite  system  of 
equations  xP  = x,  we  will  use  the  moments  of  the  stationary 
distribution.  We  presently  derive  complicated  yet  tractable 
expressions  for  these  moments. 

If  we  let  z tend  to  1 in  Equation  (7) , we  get 

oo 

(39)  X(l) (I-A)  = xn  l Bk  - x.A0. 

k=l 

Adding  X ( 1 ) n = (X(l)e)^  = (l-x0e)ii_  to  both  sides  of  Equation 
(39)  and  recognizing  that  ^(I-A+n)  = tt,  we  have 


(40) 


X (1)  = 


kL-  - 


— 1A0 


(i-A+n)  1 + (l-x  e)w. 


We  see  that  we  may  calculate  the  vector  X(l)  in  terms  of  the 
data  and  the  known  vectors  Xg  and  x-^.  This  gives  us  an 
accuracy  check  on  the  numerical  computations  of  the  addi- 
tional components  of  the  vector  x.  Having  computed  the  vec- 

k 

tors  *2'  the  sum  E ^ xy  should  be  entrywise  close 

By  evaluating  the  n-th  derivative  X^ 


to  X (1)  . 


(1)  , we 
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obtain  (1)  the  n-th  conditional  factorial  moment  of 

the  stationary  distribution  given  that  immediately  following 
a transition  the  process  is  in  state  (i, j) , for  some  i>l. 

The  quantity  X'(l)e  is  the  n-th  factorial  moment  of  the  sta- 
tionary distribution. 

Theorem  6:  The  vectors  x/1’ (1)  and  x/2^  (1)  are  given  by 


(41)  X(1)  (1)  = 


-X(l)  [l-J  kAk3  + x l B 


k=l 


+ V kB.  - x,An 
-°k=l  k -1  0 


+ (XU)  ( 1)  e)  it 


k=l 


(i-A+n) 


-1 


where 


(42) 


X (1)  ( 1) e 


2 (1-p)  l2-0 . I kBki  + 2*C>  l BkH(1)  (1) 

1 k=l  k=l 


00  hi 

+ xQ  J k(k-l)Bke  + 2x0  l kBku'  ' (1) 
k=2  k=l 


+ “°k-iBk“(2) (1)  ” 2-iAo- 


(1) 


(1) 


' XiAqu(2) (1)  + X (1) e6 (2) (1) 


- X (l)u(1)  (1)  , 


and 
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(43)  X(2)  (1)  = 


X(l)  l k(k-l)Ak  - 2X<1>(l)[l-  l kAk] 
k=2  k=l 


+ 2x  l kBk  + x l k(k-l)BkV  (i-A+n) 
k=l  _uk=2  ' 

+ (X(2)  ( 1)  e)  n , 


-1 


where 


(2) 


(44)  Xu,(l)e  = 3(i-p)  j3X(1)  (l)e6  (2)  (1) 


+ 3X(l)u(1)  (1)6(2)  (1)  + X(l)e«(3)  (1) 


+ 3x0  l k(k-l)Bke  + 6xn  l kB,_u^(l) 
k=2  K k=l 


w oo 

+ 3x0  l Bku(2) (1)  + xQ  l k(k-l) (k-2)B.e 


k=l 


k=3 


+ 3xq  l k(k-l)Bku(1) (1) 
-Uk=2  K_ 


+ 3xq  l kBku(2)(l)  + xn  ? Bku(3)(l) 
~°k=l  k~  “°k=l 


~ 3?i1A0u(2)  (1)  - xlA0u(3)  (1) 


- 1 x(1)  (l)u(1)  (1)  - X(l)u(2)  (1)  . 


Proof:  The  lengthy  derivations  are  shown  in  Appendix  II. 
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Note  that  although  the  formulas  for  the  conditional 
factorial  moments  are  complicated,  they  involve  only  known 
quantities  and  are  in  a computationally  tractable  form. 


7.  AN  ITERATIVE  METHOD  FOR  THE  COMPUTATION  OF  THE 
COMPONENTS  OF  THE  STATIONARY  VECTOR  x 

We  recall  that  the  components  of  the  vector  x 
satisfy  the  equations 

k+1 


(45)  xk  = x Bk  + l xvAk+1_v, 

\>  = 1 


for  k>l. 


We  see  that  if  the  matrix  AQ  is  nonsingular,  the  vector 
—k+1'  may  be  found  by  solving  the  appropriate  equation 

in  (45).  Neuts  has  shown  in  [27],  that  in  the  case  where  Ag 
is  singular,  the  vectors  xk+^,  k^l,  may,  in  principle,  still 
be  computed  recursively  if  the  rank  of  the  matrix  Ag(I-A^) 
is  equal  to  the  rank  of  the  matrix  [Ag(I-A^)  ^ ^ . This 
recursive  procedure,  however,  is,  except  in  very  special 
cases,  numerically  highly  unstable. 

We  suggest  computing  the  vectors  xk,  k>l,  using  the 
following  block  Gauss-Seidel  iterative  procedure. 


(46)  xk(0)  = b£ 


k-1 


Xfcln+l)  - b • + l xv (n+1) A£+1_v  ♦ ik+l(n)A0 


v = 2 


where  t>k  = (XgBk+x^Ak)  (I-A^)  and  A^  = AV(I-A^)  ^ , 
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Note  that  the  vectors  are  now  known  quantities.  Note 
also  that  in  the  recurrence  relationship  for  x^,  we  use  the 
most  recent  iterates  of  x v , v=2,...,k-l.  Using  the  moments 
that  we  have  computed  in  the  last  section,  we  first  truncate 
the  infinite  system  of  equations,  (46),  at  some  index  k* , 
where  k*  is  the  smallest  integer  not  less  than  y + 3o.  (y 
and  a being  the  mean  and  standard  deviation  of  the  queue 
length  following  departures,  respectively.)  We  continue  the 
iterations  until  the  condition 


(47) 


max 


2lkik* 


xk(n)  - xR(n-l) 


\ < 10 


-8 


is  reached.  At  this  time  we  check  to  see  if  an  adequate 
number  of  components  have  been  computed.  This  amounts  to 
computing  e,  where 


(48) 


k* 

e = 1 - I xk(n)e 


k=0 


.-4 


If  e > 10  we  increase  k*  by  1 and  continue  with  the  itera- 
tions. When  all  of  the  conditions  for  stopping  have  been 
reached,  we  utilize  an  accuracy  check  on  the  components  of 
the  vectors  xk,2*.k£k',  where  k'  is  the  number  of  components 
computed.  From  Equation  (5),  we  have 


(49)  X (1)  = l xv, 
v=l 


where  X(l)  is  known  explicitly. 


8.  APPLICATIONS 


A.  The  Mx/G/1  Queue 

We  consider  a single  server  queue  with  a general 
service  time  distribution  and  arrivals  of  random  group 
sizes,  which  occur  at  the  epochs  of  a Poisson  process.  It 
is  well-known  that  the  successive  queue  lengths  immediately 
following  departures  in  such  a queue  (denoted  by  Mx/G/1) 
form  a Markov  chain  of  type  (1)  where  the  matrices  are  all 
scalars.  In  this  case,  the  first  two  rows  are  identical 
and  the  entry  av  corresponds  to  the  probability  that  there 
are  v arrivals  during  the  service  of  one  customer.  For  the 
scalar  case, 

I=G=G  = A = 1,  and  6 = p = t 

where  A is  the  mean  arrival  rate,  £ is  the  mean  group  size 
and  a is  the  mean  service  rate.  Formula  (14)  for  p then 
simplifies  to 

1 

u = 1-p  ' 

which  is  the  classical  formula  for  the  mean  number  of  ser- 
vices during  a busy  period. 
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B.  Two  Queues  in  Series  with  Finite  Intermediate 
Waitingroom 

The  following  queueing  model  has  been  studied  by 
several  authors  [18,19,35].  A system  of  queues  consists  of 
two  units.  Customers  arrive  at  a first  unit  (I)  according 
to  a homogeneous  Poisson  process  of  rate  A.  The  service 
times  in  unit  I are  independent,  identically  distributed 
random  variables  with  common  distribution  function  H(’).  We 
also  assume  that  H(')  has  a positive  finite  mean. 

Upon  completion  of  service  in  unit  I,  all  customers 
go  on  to  a second  unit  (II)  via  a finite  waitingroom.  We 
assume  that  there  cannot  be  more  than  k customers  in  unit  II 
and  in  the  waitingroom  at  any  time.  If  upon  completion  of 
service  in  unit  I a customer  finds  the  waitingroom  full, 
then  the  unit  one  "blocks  until  a service  in  unit  II  is  com- 
pleted. At  that  time  he  is  allowed  to  enter  the  waitingroom. 

We  assume  that  the  service  times  in  unit  II  are 
independent,  identically  distributed  random  variables  with  a 
negative  exponential  distribution.  The  service  times  in 
unit  II  are  also  stochastically  independent  of  those  in  unit 
I and  of  the  arrival  process.  If  we  look  at  the  number  of 
customers  in  the  system  immediately  following  a service  in 
unit  I,  we  have  an  embedded  Markov  chain  P of  the  type  (1) 
and  state  space  {(i,j)  i>0,  1-j-k+l}  where  i is  the  number 
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of  customers  in  the  system,  who  have  not  yet  completed  ser- 
vice in  unit  I and  j is  the  number  of  customers  in  the 
system,  who  have  completed  service  in  unit  I,  but  not  yet 
in  unit  II.  The  states  for  which  j=k+l  correspond  to 
blocking . 

In  this  model,  m=n.  The  specific  form  of  the 
matrices  Av  and  Bv,  vlO  , is  complicated,  but  is  readily 
deduced  from  Formulas  (3)  - (18)  in  [18],  It  should  be 
noted  that  most  of  the  analysis  depends  only  on  the  pre- 
vailing special  structure  of  the  matrix  P and  not  on  the 
complexities  of  the  precise  definitions  of  the  matrices  Av 
and  Bv . 

C . A Single  Server  Queue  with  Versatile  Markovian  Input 

In  [30],  Neuts  defined  a general  class  of  Markovian 
point  processes,  which  generalize  the  classical  Poisson  pro- 
cess and  also  renewal  processes  of  phase  type  [23].  Such 
point  processes  are  useful  in  modelling  a large  number  of 
qualitative  features  of  arrival  processes.  Among  these  are 
group  arrivals,  randomly  fluctuating  arrival  rates  and  inhi- 
bitory phenomena. 

In  his  thesis  [34],  V.  Ramaswami  has  extended  the 
theory  of  the  simple  M/G/l  queue  to  a single  server  queue- 
ing model  having  this  versatile  Markovian  point  process  as 
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its  input  and  general  independent,  identically  distributed 
service  times. 

Although  the  detailed  definition  of  the  matrices  Av 
and  Bv  is  again  highly  involved  and  will  not  be  repeated 
here,  the  queueing  model  studied  by  Ramaswami  has  an 
embedded  Markov  chain  of  the  type  (1) . Once  the  vector  x 
has  been  evaluated,  which  may  be  done  by  using  the  tech- 
niques proposed  here,  one  can  then  draw  upon  the  detailed 
results  in  [34]  to  compute  a large  number  of  other  qualita- 
tive queue  features,  such  as  the  steady-state  distributions 
of  the  virtual  waiting  time  and  the  queue  length  at  an  arbi- 
trary point  in  time. 

D.  Queues  with  Exceptional  Services 

Consider  a queueing  situation  in  which  there  are 
occasionally  "exceptional"  services.  For  example,  the  ser- 
vice mechanism  may  occasionally  break  down,  after  which 
there  may  be  a certain  amount  of  time  needed  for  repair 
before  a service  can  be  performed.  We  can  consider  the 
breakdown  and  repair  time  combined  as  forming  an  exceptional 
service.  In  order  to  formulate  this  model,  we  shall  need 
some  of  the  basic  properties  of  phase  type  distributions 
(PH-distributions)  and  renewal  processes  of  phase  type  (PH- 
Renewal  Processes),  which  were  introduced  by  M.F.  Neuts  [22, 
23] . Only  the  basic  definitions  of  PH-distributions  will 
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be  reviewed  here.  The  interested  reader  is  referred  to  the 
cited  references  for  further  details. 


Consider  an  (m+1) -state  Markov  chain  on  the  inte- 
gers { 1 , . . . ,m,m+l } , whose  matrix  P of  stationary  transition 
probabilities  is  of  the  form 


where  T is  an  mxm  matrix  and  T°  is  a column  vector  with  m 
components.  We  shall  assume  that  the  probability  of 
absorption  into  the  state  m+1,  starting  from  any  given  ini- 
tial state,  is  equal  to  one.  This  implies  that  (I-T)  "*■ 
exists.  The  vector  of  initial  probabilities  of  the  Markov 
chain  will  be  denoted  by  and  here  we  may  assume 

that  am+l=0  * 


A probability  density  {r^}  on  the  positive  integers 
is  of  phase  type,  if  and  only  if  there  exists  a finite  sto- 
chastic matrix  P of  the  type  (1)  and  a vector  a of  initial 
probabilities,  such  that  {r^}  is  the  density  of  the  time 
till  absorption  into  the  state  m+1.  If  {r^}  is  of  phase 
type,  then  it  is  easily  seen  that 

rk  - aTk-1T°,  for  kil. 

Since  the  density  {r^}  is  determined  by  a and  T,  we  call  the 
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pair  (o^,T)  a representation  of  the  density  {r^}. 

Now  consider  the  stochastic  matrix  Q,  of  order  m, 
defined  by 

Q = T + T°A° , 

where  T?  . = T?,  for  lii,  jim,  and  A°=  diag  am)  . The 

matrix  Q is  readily  shown  to  be  the  transition  matrix  of 
the  PH-renewal  process  obtained  by  instantaneously  restart- 
ing the  chain  P after  each  absorption  by  performing  a multi- 
nomial trial  with  probabilities  a-^,  . . . ,am,  to  select  the 
new  "initial"  state.  Considering  each  absorption  as  a 
renewal,  it  is  obvious  that  the  density  of  the  times  between 
renewals  is  of  phase  type  and  is  equal  to  {r^}. 

We  can  now  construct  a model  of  a queue  with  excep- 
tional services.  Consider  a queue  with  exponential  inter- 
arrival times  with  arrival  rate  A.  We  introduce  an  under- 
lying m-state  discrete  PH-renewal  process  as  defined  above 
such  that  immediately  prior  to  a service  completion,  a 
transition  is  made  in  the  PH-renewal  process.  If  this 
transition  does  not  involve  a renewal  then  the  service  time 
of  the  next  service  has  the  distribution  F(x).  If  the 
transition  involves  a renewal,  then  the  service  time  of  the 
next  service  has  the  distribution  (x) . In  this  way,  the 
exceptional  services  correspond  to  renewals  in  the  discrete 
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PH-renewal  process.  If  we  define  Jn  to  be  the  phase  of  the 
PH-renewal  process  during  the  nth  service  and  Xn  to  be  the 
duration  of  the  nth  service  then  the  pair  (Jn,Xn)  form  a 
Markov  renewal  process  with  transition  probabilities  given 
by 

Qij(x)  = P(Jn+l=3r  Xn<x  | Jn=i} 

= TijF(x)  + TJa.jFj.U)  . 

We  define  the  matrix  Q(x)  = {Q^ j (x)  } = TF(x)  + T°A°F^(x). 

The  successive  queue  lengths  immediately  following 
departures  and  the  random  variables  Jn  form  a Markov  chain 
of  the  type  (1),  where  = Ag,  Bn  = An  for  n=0,l,...,  and 
the  matrices  An  are  defined  by 

n 

A = / e'Xu  (Xujn  dQ(x) . 
n J 

The  above  model  may  easily  be  modified  to  allow  for  a 
different  arrival  rate  during  the  exceptional  services  as 
may  be  the  case  in  certain  practical  situations. 

We  see  that  in  this  manner  a generalization  of  the 
M/G/l  arises  in  which  strings  of  ordinary  services  are  sep- 
arated by  single  exceptional  services.  The  lengths  of  the 
runs  of  ordinary  services  are  independent,  identically  dis- 
tributed random  variables,  which  may  have  an  arbitrary  dis- 
tribution of  phase  type.  Mathematically  the  queue  so 
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obtained  is  a particular  case  of  the  M/SM/1  queue.  It  is 
computationally  highly  tractable  and  permits  the  algorith- 
mic investigation  of  a number  of  control  and  optimization 
aspects,  which  we  shall  discuss  elsewhere. 

E . Bulk  Service  Queueing  Models 

1.  Bailey's  Bulk  Service  Queue.  Consider  a bulk 
queueing  model  involving  a server,  who  becomes  available  at 
the  epochs  of  a renewal  process  with  underlying  distribution 
H ( • ) • Customers  arrive  according  to  a Poisson  process  of 
rate  X.  If  k customers  are  present  when  the  server  becomes 
available,  a group  of  size  min(k,m)  enters  service.  This 
model  was  solved  by  N.T.J.  Bailey  [1]  by  the  use  of  complex 
variable  methods. 

The  successive  queue  lengths  immediately  prior  to 
the  beginnings  of  services  form  a Markov  chain  P with  the 
following  structure: 
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for  j >0 . 

We  see 

that  this 

j i 


matrix  may  be  partitioned  into  the  form  (1)  , where  Cg=A0 
and  all  of  the  matrices  Bn  and  An  are  square  matrices  of 
order  m.  We  also  note  some  interesting  consequences  of  the 
fact  that  the  first  m rows  are  identical  and  given  by  { a j } . 
We  state  these  results  in  the  following  theorem  which  may 
be  found  in  Neuts  [25]. 


Theorem  7 . For  the  Markov  chain  in  Bailey's  model,  the 
matrix  L(z)  defined  in  Formula  (18),  has  m identical  rows 
which  are  equal  to  te  first  row  of  G(z),  defined  in  Formula 
(8) . The  vector  d,  defined  in  (20)  is  given  by  the  first 
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row  of  the  matrix  G and  the  vector  d* , defined  in  (21),  is 
given  by 


d*  = w^e, 

where  is  the  first  component  of  the  vector  £.  The  vector 
Xq  is  given  by 

-1 

xQ  - v1  d, 

and 

xQe  u ^ • 

2.  Moran's  Dam  with  Infinite  Capacity.  A classical 
model,  due  to  Moran  [13],  for  a dam  in  discrete  time  with 
discretized  content,  involves  a Markov  chain  P of  the  type 
(1) , with  the  following  entries: 

c0  = Ao 
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o 

). 

= aj-i' 

if  jii/ 

= 0 , if 

j <i 

[Av 

hi 

= avm+j- 
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where  {a^}  is  the  probability  density  of  the  number  of 
units  of  water  added  to  the  dam  per  year  and  m is  the  maxi- 
mum  amount  of  water  released  at  the  end  of  each  year.  We 
assume  that  the  capacity  of  the  dam  is  infinite. 

3.  A Bulk  Service  Queue  with  a Threshold.  The  fol- 
lowing queueing  model  also  has  an  embedded  Markov  chain  of 
the  type  described  in  this  thesis.  Customers  arrive  at  a 
service  unit  according  to  a Poisson  process  of  rate  A. 

Services  occur  in  groups,  with  the  group  size  dependent  on  the 
the  queue  length  according  to  the  following  rule.  Let  there 
be  i customers  waiting  at  the  completion  of  a service.  If 
0<i<L,  the  server  remains  idle  until  the  queue  length 
reaches  L and  then  starts  serving  all  L customers.  If 
Liiim,  a group  of  size  i enters  service  and  if  i>m,  a group 

of  size  m is  served.  It  is  assumed  that  the  lengths  of  ser- 
vice of  successive  groups  are  conditionally  independent, 
given  the  group  size.  The  successive  queue  lengths  follow- 
ing departures  form  a Markov  chain  of  the  desired  structure. 
This  model  has  been  studied  by  several  authors  [5,14,17,25], 

4 . A Bulk  Service  Queue  Viewed  as  a Branching 
Process . Assume  that  a server  serves  groups  of  size  m.  If 
at  time  t=0 , there  are  i customers  we  divide  this  group  into 
groups  of  size  m,  with  any  remaining  customers  left  alone. 
Assume  that  there  are  n such  groups  of  size  m.  We  consider 
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any  arrivals  to  the  queue  during  the  service  of  any  of 
these  groups  as  the  progeny  of  that  group.  The  total 
number  of  customers  at  the  end  of  the  service  of  the  n-th 
group  will  form  the  first  generation.  If  we  continue  in 
this  manner  then  the  busy  period  starting  with  i customers 
will  be  equal  to  the  time  till  extinction  in  this  branching 
process  starting  with  i customers.  This  model  has  been  stu- 
died by  Ezhov  and  Shakhbazov  [6] . 

F.  Queues  with  Semi-Markov  Service  Times 

Queues  with  semi-Markov  service  times  have  been 
studied  by  several  authors  including  (Jinlar,  Gaver,  Loynes, 
Neuts  and  P.  Purdue  [2,7,12,15,27,29,33].  One  typical 
model  involves  an  M/G/l  queue  in  which  there  are  m types  of 
customers,  operating  under  the  first  come-first  served  dis- 
cipline. We  assume  that  the  server  expends  a random  length 
of  time  in  the  change-over  from  one  type  of  customer  to 
another.  This  model  has  an  embedded  Markov  chain  with 
state  space  { ( i , j ) , ilO,  lljlm}  where  i is  the  number  of 
customers  in  the  queue  following  a service  and  j is  the 
type  of  service  that  the  server  is  tooled  up  for  immediately 
following  a service.  etails  of  this  model  are  given  in 
[27] . Explicit  formulas  for  general  M/SM/1  queues  with 
group  arrivals  are  available  in  [29]. 
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G.  Generalized  Random  Walks 

Consider  an  infinite  Markov  chain  with  state  space 
{(0,j),  lij-n,  and  (i,j),  i-0,  l^j^m).  We  assume  that 
starting  in  state  (i,j)  i21,  l^j-m,  in  one  transition,  the 
chain  may  only  enter  the  states  { (i , j ' ) , j V j , l±j'3m,  and 
(i-l,j'),  (i+l,j'),  lij'-m}.  Starting  in  the  state  (0,j), 
only  the  states  {(0,j'),  jVj/  l-j'-nand  (l,j')f  l^j'im) 
may  be  reached.  This  model  describes  a generalized  random 
walk  and  has  a transition  probability  matrix  P of  the  form 


which  is  clearly  of  the  type  (1) . This  model  has  been  stu- 
died in  references  [33,37,39]. 

The  invariant  probability  vector  of  this  process 


may  be  evaluated  using  the  techniques  proposed  in  this 
paper.  We  also  note  that  the  matrix  P has  a structure, 
which  is  a special  case  of  a general  class  of  Markov  chains 
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studied  by  Neuts  [26]  . It  is  shown  there  that  the  invari- 
ant vector  is  of  a matrix  geometric  form. 

H.  Queues  with  Fluctuating  Input  and  Service 

Consider  a queue  where  the  arrival  process  and/or 
the  service  rate  exhibit  random  variations.  This  model  may 
be  used  to  describe  changes  in  work  shifts,  rush  hours, 
interruptions  in  the  arrival  process,  server  breakdowns, 
etc.  To  be  specific,  we  assume  that  there  is  an  underlying 
m-state  continuous  parameter  irreducible  Markov  chain  which 
governs  the  phases.  During  any  interval  spent  in  phase  i, 
the  arrivals  are  according  to  a homogeneous  Poisson  process 
of  rate  A^  and  any  service  initiated  during  such  an  interval 
has  a duration  distributed  according  to  This  model 

was  first  discussed  by  Neuts  [20].  If  we  consider  the  queue 
lengths  immediately  following  departure  we  see  that  once 
again,  we  have  an  embedded  Markov  chain  of  the  type  (1)  . 

By  making  Markovian  assumptions  on  the  service  mech- 
anism, it  is  also  possible  to  consider  models  for  which  the 
service  rate  of  a customer  may  vary  with  the  underlying 
phase  state.  This  model,  usually  called  "The  M/M/c  queue 
in  a Markovian  environment,"  has  an  extensive  literature 
[3,28,31,32,38].  It  may  also  be  treated  by  the  methods 
described  in  this  thesis,  although  results  which  identify 
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the  stationary  probability  vector  as  being  (modified) 
matrix-geometric,  provide  a more  explicit  solution,  both 
from  an  analytic  and  an  algorithmic  point  of  view. 


I . Parity  Dependent  Service  Times 

A model  which  is  related  to  the  above  "bin"  model  is 
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We  see  that  we  may  partition  this  matrix  exactly  as  in  the 
"bin"  model  to  obtain  a matrix  of  type  (1) . This  model  was 
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studied  by  Neuts  [15]  . We  note  that  the  first  row  of  the 
matrix  P2  in  [15]  should  consist  of  the  quantities  {av}, 
rather  than  { t> v } . This  also  results  in  some  obvious  changes 
in  Section  3 of  that  paper. 

J . A Bin  Model 

Consider  the  following  "bin"  model.  The  waiting 
room  of  a queueing  system  consists  of  an  unlimited  number 
of  bins  of  size  m.  As  arrivals  occur,  they  are  placed  into 
the  bins,  filling  them  one  at  a time.  When  the  server 
becomes  available,  he  serves  the  last  bin  that  received 
input,  if  it  is  not  full.  If  all  occupied  bins  are  full,  he 
chooses  one  at  random  and  serves  it.  Now  suppose  that  the 
service  times  depend  on  the  number  of  items  in  the  bin  being 
served  (i.e.  the  residue  class  of  the  queue).  The  queue 
length  immediately  following  a departure  forms  a Markov 
chain  P with  the  following  structure.  (For  simplicity,  we 
display  the  matrix  for  m=3)  . 


9.  THE  APL  PROGRAM 


The  algorithms  constructed  in  this  paper  have  been 
programmed  and  implemented  in  APL.  We  briefly  review  the 
purpose  of  each  of  the  APL  functions  .and  give  a short  dis- 
cussion of  the  more  important  ones.  In  Section  10,  we  pre- 
sent numerical  examples. 

The  index  origin  of  the  workspace  is  set  to  zero. 

The  sequences  of  matrices  {An}  and  { Bn } have  been  truncated 
at  some  index  k.  We  denote  the  matrix  A=  £ A at  the  varia- 

v = o 

ble  AA.  The  sequence  of  matrices  {Anln=0  is  represented  by 
the  three  dimensional  array  A of  dimensions  (k+l,n,m).  The 
variables  CNOT  and  BNOT  equal  the  matrices  Cq  and  Bq  respec- 
tively.  The  sequence  of  matrices  {Bn)n_^  is  the  three  dimen- 
sional array  B of  dimension  (k,n,m).  Note  that  in  the  fol- 
lowing program,  A[k;;]  = Aj.  and  B[k;;]  = Bk+^,  for  k£0. 

A . The  APL  Functions 
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Once  the  set  of  matrices  {C0 , An , Bn , n=0 , . . . , k } have  been 
initialized,  the  function  RUN  is  executed.  This  function 
calls  all  of  the  main  functions  in  the  program  and  computes 
all  of  the  desired  quantities  of  interest.  If  the  value  of 
p is  greater  than  or  equal  to  .99,  the  execution  is  halted. 
This  function  also  utilizes  all  of  the  accuracy  checks  pro- 
posed in  this  paper. 


v e e r ; a 1 i n v ; a p ; b p ;bvb  c l.  ; p ; x i j r ; r ] ; r 2 ;r  1 } 


c i :i 


1:2:1 


1:3:1 

i:  A 1 

H5::i 


1 6 1 


t:  7 :i 
1:  e j 
c?:i 

1:10:1 

fill 
11123 
L 1 3 3 
l i a :\ 


L15  1 


I"  163 
1 : 1 7 1 
1: 1 8 1 
L 1 9 3 
C203 
C21  3 
f22  3 

1:233 

C 2 A 1 


1:263 
i:279 
f 28  3 

129  1 

130  1 
1 3 I 1 


1 rs  ; e ; r 


not  , xoi-if  f ox.Jf  1 

K 1 3xSIGMAf  <X2>E  t XlE- 

A 1 I N V < fj  I - A L I f ; I 
EiVEC  r>XNOT  fl<f  0 
B P ( f B ) r A P < (f«)f  R < ■ ()  X G f ■ 1 

apl:o ; } :i<-«ro>  ;:i+.  x«iimv 

loop  1 1 ap i: « ; ; ;;i «■  a i; «<«-*  i ; ; ] + . x a i i mv 
•>  ( ( 1 If  A)  >G*  1 ) /LOOP1 
i..  a a p 2 ! e<  p i;  r b l r < r ».  j ; ; j + )Xfl  l r m v 


E ft  2 ) A 0 . 5 


M ( If  /'&  ) >R  ( 1 )/ 1.0 OP 2 


LOOP 3 J ( l<  < ;|  tf  E'P  ) / OME 
KPf  Bf  , [;  0 3 0 

APf.Ai  y to 30 

one  ; ijvl c 1',-tveo  i , ( ::oiie  i x ap|  i: <■■«  + i f ; .] ) +xhot+  , xbp[k;;;] 


: o ; ; :i 

L ) XM)  I 1 


( ( K <2  ) vK  < K 1 ) / LOOP 3 
XfH,  ( (K  : J)xM)fI 
S I A R T JIT  S f 1'  T S yl.  ; L < • 1 ..)  f p ( 1 

>'•1  v x + H > < ■■  1 1f  1 2 f MT  0 
BPECim.}  t ;|  f (Hilxh)  f lmj*  , xap 

x [;  ( m * ( l - j ) x m ) + j m | < e<  v 1 c i | ( n i ( l-  - 

> ( < k + o x * J <•  pfP+1 ) <t<  I *<  iO)fH  It 1 < ’»  If  1 2) /test 
l o o p 4 : > < I (.  I I f x [;  ( it  + ( l - jjxmj)  \ m ] { , x a p [ a ; ; j 

> ( 2 < >■><••  a • ] ) / l o o p 4 
-►SPECIAL 

r l s t ; •>  ( m a x i t x < r r s ) / o u r 
-Kl,  OOOOOOt  • U <|/|  x l x ) /ST ART 
-MO. 0001  < 1 f/x  + oxaf-K  , j )/i  nor  3 
x f • ( - M ) 4 :• 

•fO 

OUT  j 1 TOO  MAH  f ITERATIONS,  FOR  THE  VECTOR  X' 
Xf ( M ) XX 
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The  function  XVECT  computes  the  stationary  vector  x 
defined  in  Equation  (3)  by  using  the  iterative  procedure 
proposed  in  Section  7.  We  first  compute  the  index  K1  which 
is  equal  to  the  smallest  integer  greater  than  w + 3a  where 
V = X'  (l)e,  and  o = /x^2*  (l)e+X'  ( l)”e- (X'  ( 1)  e)  The  index 
K1  is  used  as  an  initial  truncation  of  the  vector  x.  We 
next  compute  the  vector  b = (bp 1 ' * * * Kl^  defined  in  Equa- 
tion (46) . The  loop  START  performs  the  iterative  procedure 
defined  in  Equation  (46).  These  iterations  are  continued 
until  Condition  (47)  is  met  or  until  a specified  maximum 
number  of  iterations  is  reached.  If  Condition  (47)  is  met, 
we  check  to  see  if  an  adequate  number  of  components  have 

been  computed.  This  is  accomplished  by  computing  the  dif- 

Kl  _o 

ference  1-^.1^  ^e.  If  this  difference  is  greater  than  10  , 

we  increase  K1  by  one  and  repeat  the  procedure. 


v Of  so  i...  v g ; «2  r ; u ; ci } t> ; d»  ; err  ; gi  ; G2  • ft 
113  ft THIS  FUNCTION  MILL  SOLVE  THE  MATRIX 
C23  (|FUHC  T I OHAL  EQUATION  G=A ( G ) 

C33  A3<-A|;»D«-D«.(  )CGITSf.03,*  n 

143  g 1 <- o x a 1 1 n v < jg  ( :t f ( t R > « , = i r < " 1 1 b ) )-A[;:L ; ; j 

C53  LOOP  1 *G<-G1  fOxGITSFGITS  + n^.Dt.tin  + OX  ( \0>f  1 1 

CA J LOOP 2 J >(2<r,  + 0x  1 1 f F 3+«2+.  XG)/1.00F2 

C73  ERR<-r/  t I G - G 1 A 1 I N V f , X A [!  0 i i 3 + «2  + ,XG+,  xG 

C83  ->(  <ci«-maxitg>gits)a<i  .000000e~8  ) <err)/loopj 

C91  ~>lt(Ci  = i)(D  ] o 15 

1103  1 A I ■ JMflXITGf  1 ITERATIONS,  THE  BEST  APPROXIMATION* 

Cl 13  1 of  G is;* 

C 123  D>°1 

C133  'WITH  A MAXIMUM  DIFFERENCE  BETWEEN  ITERATES  OF* ' (ERR 

C 1 43  >0 

C 1 53  Gf  Gl  l ((D(P,R)f  ( ( 1 - + /G1  )-r>/G2)  ) XG2FG1-G 


V 
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The  function  SOLVG  solves  the  matrix  functional 

OO  y 

equation  G = AvG  by  iterative  procedure  given  in 
Equation  (12) . i.e. 


G ( 0 ) = ( I-A-i  ) % 


(n+1)  = (I-A^  Yaq  + l AvG(n)l  . 

v=2 


The  quantity  in  square  brackets  is  evaluated  by  Horner's 
algorithm  for  the  computation  of  a polynomial.  As  is  well- 
known,  this  algorithm  minimizes  both  the  number  of  matrix 
multiplications  and  the  number  of  matrix  summations  required. 
We  continue  the  iteration  until  a specified  maximum  number 
of  iterations  have  been  reached  or  until  the  condition 


lii, j-m  | Gi j (n+1) -Gij (n>  | < 10  ' 


has  been  reached.  Finally  a linear  extrapolation  is 
applied,  which  sets  the  final  matrix  G equal  to  (G^j),  where 

Gij  = Gij<n+D  + 9i  [Gi  j (n+1)  -G^j  (n)  ] , for  l5i,jim. 


The  quantities  0^  are  determined  so  that  the  row 
sums  of  G are  equal  to  one. 


V F.T<-STVECT  PfMJPMjMJ 

Cl]  n THIS  FUNCTION  WILL  CALCULATE  THE  ST  AT I OMART 
C23  flPROBflBILm  VECTOR  OF  THE  STOCHASTIC  MATRIX  P 

C33  P«~<  <PfFM)  ft'  >PMK‘l»-l  + HMf/>P  )fP  )fF 

C43  PIpPM+,X0(  ( HfPM ) + < lMl)o,  = lMlfltMFfP)-F 

C5II  PIf (PI f 1-+/PI 
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The  function  STVECT  calculates  the  stationary  pro- 
bability vector  of  an  irreducible  stochastic  matrix  P,  using 
using  a method  suggested  by  P.  Wachter  [36],  We  write  the 
stationary  equations  _nP  = _n,  as 
m 

j ( 6 ^ j — P j i ) = 0/  for  i=l,...,m. 

where  is  the  Kroneder  delta.  If  we  add  Pm^  to  both 

sides  of  the  i-th  equation,  we  have 

m 

1 ’j<sij-pjitpmi>  * pmi- 

3=1 

Note  that  the  first  m-1  equations  do  not  involve  the  quan- 
tity irm.  Wachter  has  shown  that  the  first  m-1  equations, 

m-1 

1 n j ^ '"’i  j-p  ji+pmi^  = pmi» 

]=1 

form  a nonsingular  system. 

In  matrix  notation,  we  define  I to  be  the  identity 
matrix  of  order  m-1,  the  vector  P^  to  be  the  (m-1) -vector 
whose  i-th  component  is  Pmi , the  matrix  P*  to  be  the  matrix 
obtained  by  deleting  the  last  row  and  last  column  of  P,  and 
the  matrix  Pm  to  be  the  matrix  all  of  whose  rows  are  identi- 
cal and  equal  to  P^.  The  above  system  of  equations  may  be 
written  as 

n*(I-P*+Pm)  = P,n, 
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where  jr_*  = ( irj_ , n2  , . . . , nm_j_)  • 

The  function  STVECT  solves  this  system  for  and 
then  computes  nm  by 

m-1 


V MUfMUVECTfS 

till  0 THIS  FUNCTION  CALCULATES  THE  VECTOR  MU 
C2!  CiTILDAf  ( f G ) f'QVECT (-STVECT  G 

C31  Zf(DELTABfDIAG  BETA)+, xGTILDA 

C4  J MU<-  f/  ( GTILLA+I-G  ) + , X0(  I-AA  ) +GTILDA-2 

V 

The  function  MUVECT  computes  the  vector 
u=  (I-G+G)  { I-A+G-A  { g ) G) --Le . 


V 

k f k matr  i ; r ; i e< 

til 

( ) R ) o , = | Rf  Iff  BNOT 

C21 

i 

V 

K(-APRIME+(CPRIMEf(CNOT+,  X0IB-BNOT)  ) + , xBFRIME 

V 

Lf LMATRIX 

Cll 

< AAPRIME(-0I-APRIMEf  ASERIE S ) + , xCNOT 

C21 

Lf  BHOT-f  ( BPRTMEfBSERIES)  + , XH 

y 

The  functions  KMATRIX  and  LMATRIX  compute  the 


matrices  K and  L respectively,  using  their  defining  Equa- 


tion 

(19)  . 

V 

HSTARf SOLVH f Z f C ( T 

Cll 

z«-(  1+f «)f  0xC<-2 

C21 

LOOP  ; -> ( <~i  + ifpA)  >c<-c+(  \0)f  ltl , 

F Z«-Z+CxA[C+l  f ; ])/LOOP 

C3J 

ZfZ+A[2f t J 

C41 

T<-  ( ( AA-A[  Of  f J )-APRIME-AAAfZ+  , 

XGTILDA ) + , X INVERSMU 

C5.1 

HSTAF:pAAPRIME+  , X 1+T 

v 
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Cl] 
C2] 
C 3] 


KSTflKfSOLVK  ( Z J E<E<E< 

BEHM-BTILDAo  , xOVECT 

Z<-  ( ( ftA~ft[;o;  ; ] ) - APR  I ME  ) +AAA+CPRIME+  , X ((+/*■  ) -BPRIME-6BB  ) 
K5TARfl+ ( +/CPKIME )+2+ , X IHVERSMU 


V 

V 


Cl] 

C2] 

C3] 


DSTARlfSOLVL ( Z ) BGTIUDA 

BGTILBAf  ( GMUf-GVECT+  , X MU  ) X BT  I LDOf  ( +/B1MEAH  ) -1-  + /BHOT 
Zf.  ( ( f /?!  ) + , X IHVERSMU*-  ( 0GTILDA  + I-G  ) + , XMU  ) +EGTILHO 
tiSTflRlH+(BPPIME+,X  (HSTARfSOLVH) -IHVERSMU  ) +Z 


The  functions  SOLVH,  SOLVK,  and  SOLVL  compute  the 
vector  h* , <* , and  d* , respectively.  These  vectors  are 
defined  in  Equation  (22) . 


V 

Cl] 

<7 

V 

Cl] 

C2] 

V 

r<2<-  r>  e i_  t a 2 

D2fPI+,  X ( + /A2<“«2MErtH  ) + ( 2x  A1  + . XUl  ) -RHOxUJ 


D34-DEUTA3  ; Z 

Z<-  ( +/A3*.A3MEAM  ) - ( U 2 X R H O ) f E'  2 X U 1 
B3ffi  + .x  ( 3XA1  + , XU2)K3X«2+.  XUl  )+z 


The  functions  DELTA2  and  DELTA3  compute  the  quanti- 
ties 6^2Nl)  and  6^(1),  respectively,  using  the  recurrence 
relations  given  in  Equation  (31). 


V U 1 f-U  1 F'R  I M 

Cl]  lJl<“  ( ( «1  I H V *■  Q I — A A — ( f A A ) fPI  J + .XBETA)  -RHC) 

V 


Cl] 

C2] 

C3] 


v u 2 <-  u 2 F ' R i m j k i f k 2 

K l«-2x«l  IHV+  , X ( A 1 —RHC)  X I ) + , X U 1 *-U  1 PR  I M 
K2f  4 /A1  IHV  f , x ( A2-  ( D2*  tiEUTA2  ) X I ) 

U2<"K  1 + K2-2X  (VlfVlPRIM)  + , XUl 
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V u3fU3FKiM ; K 1 ; z 

LI  1 1 2*-  ( -3x»2xui ) + ( +/«3 ) -D3«-r.ELT«3 

[23  k if«i IHV+ , x(3x«l+. xU2)+(-3xKHOxU2)+(3xA2+. xui )+z 
13  ] u3«-ki-3x  (vi  + ,xU2)  + (V2f-V2FRiM)  + .xui 

V 


v VlfV]PRIM 

c 1 3 V1(-(FI+,X»1+,X«1I  HV  ) -RHO  X P I 

V 


v v2<-v2F:ri::iM ; s 
Cl  3 2+-2xf:hoxvi<-viprim 

C23  V2*< (F i+, x«2)+(2xvi+, x«l )-(D2xpi )+2)+, xAiiMV 

v 


The  above  functions  compute  the  successive-  deriva- 
tives of  the  Perron-Frobenius  eigenvectors,  u^^  (1)/ 


(2) 


(1),  u ( 3)  (1)  , vUJ  (1)  and  v ^ 2 ^ (1)  , respectively.  These 


,(1) 


(2) 


derivatives  are  also  defined  in  the  recurrence  relations 
given  in  Equations  (30)  and  (31). 


i 


V M]E(-K]MEANJ  K ] £ K 2 J 2 

Cl  3 «/]<-  < XHOT+  , X+/R  ) -XONE+  , X«C0  » T 3 

C23  '<!<-(  < 1-+/XNOT)  xt'2 ) +«4+  * x ( ( 2 x u i > +u2«-U2PRIm  ) 

C33  2f(+/xNOT+,  XB2<-B2ME«M)  + K1 

C 4 3 k2K2x::mot+,  x ( bi«-e>1meam ) + . x ( 1+ui  ) ) +z 

C53  X1U1F>(XK-(«4+.XAHHV)  + (1-  + /XWOT)xFI)  + (xU1 

C63  K2-r2xl-RHO)-xiuiP 

V 


C 1 3 
C23 


XJ  P<-K;lF  RIM  ; K 1 
K 1 «-  ( X1E«.  XJMEON  ) xFI 

X 1 F <-  K ] *•  ( ( «4  + XHOT+  , XB1  )-!<!+,  X ( I-«l  ) ) + , XAIIMV 


v X2E«-:''2MEOWf  Kt  > K2»  K3$y$z 

C13  2<-<  (+/»)  + . x ( K2«-(U3<-U3PRiM)+3xU2)  ) -m-/»3<-*3me«n 

C23  KifXMOT+ , x (3x»2+ ,x(i+ui> )+(3xfl+. x(u2+2xui ) )+z 
C 3 3 ■»  < 3 x t'2  x >•:  1 + . x u i ) + < i _+/xnot  ) x r>3 

C43  K3K-XOME+,  x«C0?  i 3+  . x K2 ) + ( 3xE'2xS<lE ) + »' 

C33  >:2E«-<  ( « l + «3 ) + 3x  1 -Rho  ) - < xi+  , xU2 ) +4x«lF+ . xui^-3 

v 


A 


I 
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v 

C1J 
C23 
C31 

v 

The  functions  X1PRIM  and  X2PRIM  compute  the  moment 

vectors  X' (1)  and  Xv  (1),  respectively.  The  functions 

X1MEAN  and  X2MEAN  compute  the  quantities  X' (l)e  and 

(2) 

X (l)e,  respectively,  which  are  needed  in  X1PRIM  and 
X2PRIM.  All  of  the  above  quantities  are  defined  in 
Theorem  6. 


v flF  RIMEf  flSEP:IES  JC 
LI]  flPRiMEFft[CK (f«)-l)C0D5 n 

L2]  LOOF' }-»(2<C  + 0X  1 1 |OFRIMEffl[CfC-l  ; ( ]+flPRIME+  , xG)/LOOP 

L3]  OFF:IME<-AL  1 i f ]+AF  F:IME  + , XG 

V 

V PPP:IMEfS5ERIE5;C 

LI]  BFRIME«-BLC«-(  (fB)-l)LO];?l 

L2]  LOOP*  + < 1 <C  + 0*  1 1 fBPPIMEfBLCfC-1  ; } ] + BPRIME+,  XO)/LOOP 

L3]  EPP.XMEpBLO  f 5 ]+BPRI  ME+  , XG 

V 

The  functions  ASERIES  and  BSERIES  compute  the  quan- 
v V_1  r v-1 

titles  l AyG  and  l BVG  , respectively,  using  Horner  s 
v=l  v=l 

method . 

V fllfAlMEAiqC 

Cl]  «:L<  < 14-f «>fOxc«-i 

L23  uooF  j«i<-ai+alc;  ; ]Xc 

C33  •♦(  (Iff  A)  >c<-c+d/loof 

V 

v a2<-A2MEam;c 

«2«-(  14-f  «)f  OxCf-2 
loop  • A2«-A2+alc;  ; ]Xcxc-i 
•♦<  < Iff  «)  >C<-C+1  ) /LOOF* 


::2F'<-:':2FFi:IM ; k i ; z 

K 1<-  ( K2E<-K2MEAt(  ) XPI 
2PXHOT+,  XB2  + 2XE-1 

X2p<-«1  (X1  + , xA2)+s-2x::if  + , x ( i-«l ) ) + , x«1ihv 


HID 

T2J 

C33 


v 
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v «3<-A3MEr«w;c 
T1D  «3<-  < 14-f  a)>'0xc<-3 

C23  l.oof  :a3<-a3+acc;  ; ]Xcx  <c-i ) xc-2 
C33  -»  ( ( 1 ) > Cf  C.fl  ) /loop 

v 

V r<  1 f • E<  1 MEAN  ; c 

C 1 3 E<i<- ( 14-P»)pc'«-o 

[23  LOOP  ‘Blf.Bl+B[C(  f 3XC+1 

133  ■»(  (it/' 6)  >Cf.C+l  )/LOOF 

V 


v E<2<-»2MEftW i c 
C13  *2<-(  14-f  *)f  Oxc<-i 

[23  LOOP  }E<2Pf2+E,Cprf3xcxCfi 

C33  ■»<  ( lfpB)  >CfC  + l )/LOOP 

V 


v &3pE‘3meamjc 
C 1 3 E‘3«-(14-fE‘)poxCf.2 

[23  Loop;£>3fB3+B[C}  » 3xcx  < c + l )xc+2 
C33  ->(  (Iff »)  >c<-c+  d/loop 

v 


The  functions  A1MEAN,  A2MEAN,  A3MEAN , BIMEAN, 


B2MEAN  and  B3MEAN  respectively  compute  the  moment  matrices 


k 

l vV 

v = i 


k 

y vb  , 

£,  v 
v=l 


v = 2 


v (v-1) Ay , 


v = 2 


v (v-1) Bv , 


k 

y v(v-l)  (v-2)Av, 
v = 3 


k 

l v(v-l) (v-2)Bv. 
v=3 


V DEL.TOPf.EH  AG  PJN 

C13  n THIS  FUHCTIOH  WILL  CREATE  A DIAGONAL  MATRIX 
C23  n WITH  THE  ELEMENTS  OP  P ALONG  THE  DIAGONAL 
l 3 3 DELTAPf  Nf  t<  , ( Nf-  ( f P ) f f P ) p 0 

V 
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The  function  DIAG  creates  the  matrix  Me)  which  is 
a diagonal  matrix  with  the  elements  of  the  vector  £ along 
the  diagonal. 


B.  The 

Global  Variables 

A 

The  following  is  a list  of 

- A three  dimensional  array 

A1 

n — 

k 

~ I.  vAv 

A2 

V = 1 

k 

- 1 v ( v-1) Av 

A3 

v=2 

k 

- 7 V (v-1)  ( v-2) A 

A4 

v=3 

k 

“ *0  l Bn  - *iA0 

A1INV 

n=l 

- (i-A+n)-1 

AA 

k 

- A = l A 

AAA 

v=0 

k 

- 1 (v-l)AvG 

AAPRIME 

v=2 

k , -I 

" (I  - l A VGV_1) 

APRIME 

V = 1 

k i 

- ? 

V = 1 

An  for 
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B 

A three  dimensional 
N=0 , . . . , k-1 . 

Bl 

- 

k 

l vBv 

V = 1 

B2 

- 

k 

I v(v-l)Bv 
v = 2 

B3 

- 

k 

l v(v-l)  (v-2)Bv 
v = 3 

BBB 

- 

k - * 

l (v-1) BmG  = 1 vB 

v=2  v=l 

BETA 

- 

6 

BNOT 

- 

B0 

BPRIME 

- 

r v-1 

I B^G 

V = 1 

CNOT 

- 

n 

o 

CPRIME 

- 

c0(i-Bo)_1 

D 

- 

d 

D2 

- 

6 (2) (1) 

D3 

- 

6 (3) (1) 

DELTAB 

- 

A (0) 

DSTAR 

_ 

d* 

BN+1' 
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G 

GITS 

GMU 

GTILDA 

GVECT 

H 

HSTAR 

I 

INVERSMU 

K 

KAPPA 

KSTAR 

L 

M 

MAXITG 


- G(l) 

- The  number  of  iterations  required  to  compute  the 
matrix  G. 

- 2iL 

- G 

- a 

- H ( 1 ) 

- h* 

- An  m*m 

- (I-G+G) 

- K ( 1 ) 

- 

- <* 

- L (1) 

- The  order  of  the  A-matrices. 

- The  maximum  number  of  iterations  allowed  for  the 


identity  matrix 
-1 

P 


computation  of  the  matrix  G. 
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MAXITX  - The  maximum  number  of  iterations  allowed  for  the 
computation  of  the  vector  x. 

MU  - jj 

N - The  order  of  the  matrix  Bg. 

PI  - £ 

RHO  - p 

SIGMA  - a = /x"  (1)  e + X'  (1)  e - (X1  (l)e)2 

TIME  - The  CPU  time  used  in  executing  the  program. 

U1  - u(1)  (1) 

U2  - u(2)  (1) 

U3  - u(3)  (1) 

VI  - v{1)  (1) 

V2  - v(2)  (1) 

X - -x 

XI  - X(l) 

XlE  - X’ (1) e 

X1P  - X' (1) 
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X1U1P 

- X(l)u(1) 

X2E 

- X(2)  ( 1)  e 

X2P 

- X(2)  (1) 

XNOT 

" -0 

XONE 

" *1 

C.  The  Structure  of  the  APL  Program 


We  present  below  a tree  diagram  representing  the 
structure  of  the  APL  program.  We  display  the  order  in  which 
the  functions  are  called  as  well  as  where  each  global  varia- 
ble is  first  initialized.  A rectangular  box  represents  the 
function  being  called,  an  oval  represents  the  resultant  of 
that  function  and  a diamond  represents  any  global  variables 
which  are  initialized  internal  to  the  function.  We  assume 
that  the  variables  A,  AA,  B,  BNOT , CNOT,  M,  N,  MAXITG,  and 
MAXITX  have  already  been  initialized. 


RUN 
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10.  SOME  NUMERICAL  EXAMPLES 

We  demonstrate  the  use  of  the  algorithm  discussed 
in  this  paper  by  some  numerical  examples.  As  we  noted  in 
Section  8,  in  many  applications,  the  detailed  structure  of 
the  sequences  of  matrices  {An}  and  (Bn)  may  be  quite 
involved.  In  order  to  avoid  cumbersome  numerical  integra- 
tions, we  chose  to  present  only  simple  examples.  As  such, 
their  practical  significance  may  be  somewhat  limited,  but 
they  do  illustrate  the  concepts  discussed  in  this  paper. 

Consider  a queue  with  deterministic  service  times, 
i in  which  a service  time  of  length  c^  is  followed  by  m-1  ser- 

vices of  length  c.  The  next  service  is  again  of  length  Cj_, 
and  the  pattern  is  repeated  periodically.  During  a service 

of  length  c^,  there  are  Poisson  arrivals  of  rate  X.  and 

**■ 

during  services  of  length  c,  there  are  Poisson  arrivals  of 
rate  X.  It  is  easily  seen  that  in  this  case,  Cg  = Aq , and 
that  the  sequence  of  matrices  {Bn}  is  equal  to  the  sequence 
of  matrices  {An>.  The  mxm  matrices  An  are  given  by 
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An  = 


-X  c 

e 1 1(X1c1)n 


e~Xc(Xc)n 
n ! 


n 


e~Xc(Xc)n 


— Xc 

e Uc)n 
n! 


. . 0 


. . 0 


Since  A = £ An  is  doubly  stochastic,  its  invariant  proba- 

n=0 

bility  vector  £ is  given  by  n = — e'.  The  vector  8 is  equal 

- ■ m ~ — 

to  (X^c^,  Xc , . ..,  Xc) , and  therefore 


p = tt 3 = -[X.c,  + (m-1)  Xcl  = Xc  + 
— m ii 


Xlcl  " Xc 
m 


In  the  numerical  examples,  we  further  simplify  by 
assuming  that  the  arrival  rate  is  constant  during  all  ser- 
vices, i.e.  X j_=X , and  normalize  X to  be  equal  to  one.  The 
matrices  An  are  generated  by  the  following  program: 
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v a « ( j f.  ■ 1 1 k:  r n r e ; » e r ; m ; z 

til  i c 1 <•  i.  a m r>  ft  i x c j 

I' 2 ] LCMflMDftxC 

E3J  BETf-LCl  , ( M — 1 )fLC 

L 43  2<-0X«<-J  etilftG  Mf  (*~L.C)xM<-l 

L. 5 1 «[  Of  1 ]<•  < * -LC1  ) 

C61  LOOF  :«(  n,[--l  + HiH‘t]lelUflG  Mf  ( *-LC  ) x ( LC«N  ) ~ { M 
r 7 1 Of l ]<  ( *-LCl  ) X(LCl*H)iJ N 

coi  jf2+Hxnc«*n 

L?:i  >(  (>lt  )„<4)/L.OOF 

CIO]  -><  ( r /(*='£■ ' - (+/?)  «Nxl V/I/A)  ) > 1 « 000000E”8  ) /loop 

C 1 :l  .1  « f « t C 0 3 ( A « <•  1 & " 1 A ° h f 1 ) - + / « 

C 1 J V.-!U)1  f-CHo  r f-^rof  f ] 

ct.n  i o o 

V 


Assuming  the  parameters  LAMDA,  LAMDAl , C,  and  Cl  have  been 
initialized,  the  matrices  An  are  generated  until  the  maximum 

n n 

element  in  the  vector  £ - £ vA  e - (n+l)f~e  £ A e'lis 

v=0  v L v=0  V-J 

— 8 

less  than  10  . That  this  is  an  adequate  truncation  has 

been  shown  in  Neuts  [24].  The  last  matrix  is  added  to  the 
sequence  to  ensure  that  the  computed  matrix  A is  stochastic. 
We  ran  the  program  several  times  with  m=5,  varying  the  para- 
meters C and  Cl. 


The  results  of  the  first  run  are  shown  with  a com- 
plete listing  of  the  output  and  the  remaining  runs  are 
shown  with  an  abbreviated  listing.  Also,  to  conserve  space, 
we  only  print  5 digits  although  AFL  computes  to  18  digits 
and  all  of  the  accuracy  checks  hold  to  13  digits  at  least. 
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COMPUTATION  OF  1 IT  CL'  STATIONARY  DISTRIBUTION  OF 
THE  QUEUE  LENGTH 

A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A 


M E 

QUALS  5 

THE 

ARRIVAL. 

RATE 

LAM  DA 

THE 

ARRIVAL. 

RftTE 

L..AMx:*n } - 

T ME 

SERVICE 

1 X ME 

C-  =•■  0.1. 

THE 

SERVICE 

TIME 

c-l  = 1.6 

THE 

MATRIX  A 

IS  • 

0 1 0 0 0 

0 0 .1  0 0 

0 0 0 1 0 

0 0 0 0 1 

1 0 0 0 0 

THE  ROW  SUMS  OF  A ARE ♦ 

1 :i  1 1 l 

A A A A A A A A A A A A A A A A A A ft  A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A 

C O M P U r A T I O N O F T H E AU  X 1 1...  L A R Y « U A M T I T I E S 

A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A 

THE  VECTOR  PI  IS  J 

0.2  0.2  0.2  0.2  0.2 

T H E V EOT  O R B E 7 A . I S J 

1.6  0.1  0.1  0.1  0 . 1 


RHO  EQUALS  0.-1 
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i mp:  matrix  r i s * 

0.1540? 
O.oj 1177 
0 . 0022:^ 
0.90701 
0.01 52 1 8 

9 ITERATIONS  WERE  EE  CM  1 1 RE V FOR  THE  COMPUTATION 
O F THE  MATRIX  G 

THE  VECTOR  G ISJ 

0.20625  0.065771  0.15579  0.22574  0.26645 

THE  VECTOR  MU  ISJ 

2.074  1,1118  1.115  1 . 1747  1.2797 


0.078351 
0.0010222 
0.011303 
0.0030 15 
0.91275 


0 , 20972 
0.000074:1  1 
0.00041068 
0.0025971 
0.019074 


0.3061 
0 . 90495 
0.00064713 
0.0039732 
0 . 020702 


0.25174 
0.081974 
0 . 9054 
0 , 0034092 
0 . 024258 


ft  ft  ft  ft  ft  ft  ft  ft  A * ft  ft  ft  ft  A ft  ft  ft  ft  ft  ft  ft  ft  ft  A ft  A ft  ft  ft  ft  ft  A A ft  ft  ft  ft:  * ft  ft  ft  A ft  ft  ft  ft  A ft 
AAAAAAAAAAAAAAAAAAAAAAAAAAftAAAAAAAAAAAAAAAAAAftAAA 

AH  ACCURACY  CHECK  ON  THE  COMPUTATION  OF  THE 
MATE  IK  © AH  1.1  THE  VECTOR  MU  AS  SHOWN  IN  COROLLARY 
l IN  THE  THESIS  IS  THE  FOLLOWIHG* 

THE  INNER  PRODUCT  OF  THE  VECTORS  © AMD  MU  IS* 

1 .66  6 6 6 6 6 6 5 5 

THE  CM  ) A H T I T Y ->  % ~ ( ;|  _ R H O ) I S * 

1 . 6 6 6 6 6 6 6 6 5 5 

ft  ft  ft  ft  ft  A ft  A ft  A ft  ft  A ft  A ft  ft  A ft  ft  ft  A ft  ft  A ft  ft  ft  ft  A ft  A ft  ft  ft  ft  ft  ft  ft  A ft  ft  ft  ft  ft  A A ft  A A 
A ftftftftftftftftftftftftftftft.ftftftftftft.  ftftftftftftftftftftftftftftftftftftftftftft  ftftftftftft 


THE  VECTOR  t<  IS* 

0.28625  0.065771  0.15579  0.22574  0,26645 

THE  VECTOR  KAPPA  ISJ 

0 . 1. 1 0 3 0 . 2 9 7 0 . 2 6 1 5 2 0 . 1 9 1 8 8 0 , 1 3 9 3 
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I 


' 

i 


r h e:  v e:  c r o r d s r a r i s ; 

2 . 8 7 4 1.1 1 1 8 1 . 1 1 5 1 . 1 3 A 7 1 . 2 7 V 7 

T HE  V E C T O R H S T A R I S ‘ 

2 . 874  1.11 18  1 . 1 15  1 . 1347  1 , 2797 
the:  vector  kappa-star  is; 

3.0361  5.5253  4.9953  4.4094  3.7609 

************  ft  A**  ft  ***********  A ************  ********* 

COMPUTATION  OF  THE  INVARIANT  VECTORS  KNOT  AND  NONE 

************************************************** 
TH E V E C T O R N N O T I S • 

0. 17175  0.0 3 9 4 6 3 0 . 0 9 3 473  0.1 3 5 4 4 0 . 1 5 9 8 7 

THE  VECTOR  NONE  IS* 

0.023709  0.06384  0.056213  0.041245  0.029943 

************************************************** 
************  A ***********************  A ************  * 

THE  ACCURACY  CHECK  PROPOSED  IN  COROL.L.AR  r A 
VERIFIES  THE  COMPUTATION  OF  A EL.  I HE  QUANTITIES 
INVOLVED  IN  COMPUTING  THE  VECTORS  NOT  AND  NONE, 

THE  VECTOR  KNOT  SHOULD  BE  EQUAL,  TO  THE  QUANTITY' 

KNOT  TIMES  BOOT  PLUS  NONE  TIMES  CNOT 

THE  ABOVE  QUANTITY  ISJ 

0 .171 75  0 . 039463  0 . 093473  0 . .13544  0 . 15987 

*************************************************** 
*****************  A A **********  A A * * A A AAA  A A A A A * A A A A A A A 
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COMPUTATION  OF  THE  MOMENT!-;  OF  THE  S T A 7 I OMARI 
VECTOR  X 


ft.  ft  ft  ft  A A * ft  ft  ft  A ft  ft  ft  ft  ft  ft  ft  ft  A ft  A A * A A ft  ft  ft  ft  ft  ft  A ft  A ft  A A ft  ft  ft  ft  ft  ft  A A A ft  A ft  ft 


THE  CONDITIONAL.  MEAN  HUEUE  LENGTHS  IN  THE 
V A R I O U S F H A S E S A R E ♦ 


0.17 1 04  1 . 6 2 9 8 0 . 9 2 7 1 1 0.4944 8 0 . 2 7 1 6 8 

THE  MEAN  (QUEUE!  LENGTH  I S J 


0 . 69882 


THE  STANDARD  DEVIATION  OF  THE  QUEUE  LENGTH  IS* 

1 . 2825 


ft  ft  ft  A ft  ft  ft  ft  ft  ft  ft  ft  ft  A ft  ft  ft  ft  A ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  A ft  ft  A A A A ft  A A ft  ft  A A A ft  ft  ft  ft  ft  A 
THE  STATIONARY  VECTOR  X 

ft  A ft  A ft  A A ft  A ft  A ft  A ft  ft  ft  ft  A A A ft  ft  ft  ft  ft  ft  ft  ft  ft  A ft  A ft  ft  ft  ft  A A A A ft  A ft  ft  ft  A A A ft  A 


SO  COMPONENTS  OF  THE  VECTOR  HAVE  KEEN  COMPUTED 
AND  THESE  COMPONENTS  ADD  UP  TO  0.99996 

the:  vector  x in  partitioned  form  is  given  by; 


0 

1 . 7175E“1 

3 , 9 4 6 3 ® " 

2 

9 , 34 73® “2 

1.3544®“! 

1 .5987®“i 

1 

2 . 2 7 0 9 e ""  2 

6. 3840 ®“ 

2 

5. 62 13® “2 

4. 1245® "2 

2. 994 3® ""2 

n 

3 . 4 6 53® “3 

5. 1794*' 

3, 06 13®  "2 

1 .5423®  "2 

7.2202®  "3 

3 

B.  04  5 9'  ®“  4 

2.8137®' 

o 

1 . 3235® “2 

5 . 5538®' “3 

2 . 158 7® “3 

4 

2. 04  24  ® “4 

1 . 1537E" 

P 

4, 6364® “3 

.1. . 7238® '"3 

6. 0561®'"  4 

5 

4 . 8375® ' 5 

3 . 8205®"' 

3 

1 .3 695® “3 

4.6585®  4 

1 .5238® “4 

6 

9 .09 49®'“  6 

.1 . 0689*" 

3 

3 . 5236®'"  4 

1 , 1.219®  "'4 

3 . 4791®  ""5 

7 

1 , 2 52 4 e~ 6 

2.6139®"* 

4 

8. 1079® “5 

2.4590®'  5 

5 . 6996®' “6 

8 

1 ,21 35® “7 

5,7436®“ 

5 

1 . 7057®  " 5 

3. 1938® “6 

6 . 1207® "7 

9 

7 . 2033®”9 

1 . 1606®“ 

5 

1 .3537® “6 

2. 1.333®"  "7 

3 . 7920® "8 
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****************************************  ft  A ********* 
*************************************************** 

THE  ACCURACY  CHECK  PROPOSED  If!  THE  THESIS 
F 0 1. L O W I N G E « IJ  A T I O N (40)  I ->  V E R I F I E D A S 
1-01.1.  OWS  ♦ 

THE  MAXIMUM  DIFFERENCE  BETWEEN  THE  VECTOR 
GENERATING  FUNCTION  X EVALUATED  AT  ONE  AND 
THE  SUM  OF  THE  PARTITIONED  COMPONENTS  OF  X 
( E X C L..  U D I N G X N O T ) IS* 

e.8925E"6 

***************************  A A A * ******  A ************  ft 
*************************************************** 
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computation  op  the  stationary  DisTRjeurioM  op 
THE  QUEUE  LEHGIH 

A************************************************* 


M EQUALS  5 


THE 

ARRIVAL 

Ei:  A T e;: 

LAM  DA 

- 1 

THE 

ARRIVAL 

R A T Ei! 

LAMD A 1 

“ 

T PIE 

SERVICE*: 

TIME:. 

C s o. 

49 

THE 

SE  Ft  VICE 

T I M lii! 

(•■1  " 0 

.04 

RHO  EQUALS  ()  » 4 

THE  VECTOR  mu  is  1,0754  1.8067  1.0458  1.7 70S  1.6078 

******************  Art  A 

T H E V Ei:  C T O Fi:  X M O T 1 S ♦ 

0.10467  0.16272  0.11698  0,10937  0.10626 

T H E V E C T O R N 0 1 1 E I S ♦ 

0.06469  0.028 2 2 2 0.0 6 1 5 4 5 0 . 0 6 4 0 9 0 . 0 6 4 5 9 5 

ftftftftftftftftftftftftft ft ftftftftftftftft A **************  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft 

THE  COME*  I TJOMAL  MEAN  QUEUE  I.EMGTH5  IP)  THE 
V A R I O U S F'-  PI  A S E S A P:  E • 

0,68219  0.24554  0.54915  0.62403  0.66086 

THE  MEAN  QUEUE  LENGTH  IS* 

0 . 55235 

the:  STANDARD  DEVIATION  OF  THE  QUEUE  LENGTH  IS, 


0.84426 
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ft  A A A A A A A A A'  A A A A A A ft  ft  A ft  A ft  A A ft  ft  A A ft  A A A ft  A A A A A ft  A A A ft  A ft  ft  ft  A A ft  ft 

THE  VECTOR  K XM  PAK  T ITIOMEt'  FORM  13  GIVEM  T r ; 

1.  , 1A9CJE “1  1 . 0937E“"  1 .1. . 0A2AE  "'1 
A . I 5 A 5 £::  2 A . 4 0 9 0 e 2 A . 4 5 9 3 £ ' 2 

I . • I.  38E ■’  2.044OE-  2 2. 1873E'"2 

3 . 5 7, 5 0 e - 3 4 . 9 7, 4 3 e “ 3 5 . 6 5 7 7 e " 3 

6 . 4993E  "•  4 9 , 966!tH  4 1 . 2741. E“ 3 

1 . .1  9 2 0 e - 4 1 , 0 9 9 4 e 4 2 . A 5 1 1 e ~ 4 
2 . 1 1 0 9 e - 3 . 4 7 0 5 e ~ 5 5 .205 5 e ' 5 
3 . 29  4VE  ,<>  5 . 78B7E-6  9 . 1 23AE~A 


A ft  ft  ft  ft  ft  A ft  A A ft  A ft  ft  ft  ft  » ft  ft  i ft  ft  A A ft  ft  ft  ft  ft  ft  A ft  ft  ft  ft  A ft  ft  ft  ft  ft  ft  ft  ft  A ft  A A A ft  ft 


0 1 . 0467E  1 1.  , A272E  1. 

1 6 . 4 A 9 0 e ” 2 2 , 0 2 2 2 e ~ 2 

2 2 . 2S98E"2  A . 09 79 e ' 3 

3 A . 1 400E  ""3  1.6724E”3 

4 1.4 75 I E ~3  3 . 786 9E  " 4 

5 3.299AE -4  u . 0AA4E  5 

6 A.  951  OK 5 1 . 5278E'"  5 

7 1.2 8 4 1 E -•  5 5 . 5 0 4 A e “ 7 
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COMPUTATION  OF  THE  STATIONARY  DISTRIBUTION  OF 
THE  CUIEUE  LENGTH 

* **ft  #»*****»  M * * * * A A H * A *******  )l  ************  A * H * 


m el 

GUAL..S  5 

the;: 

A RftIVfiL 

RATE 

LAM  DA  :::: 

the 

APR X VAL 

R A T EL 

LAML'Al 

T HE 

SERV I CE 

t :i:  me 

C ~ 0.4 

T HE 

SERVICE 

Y x m e: 

< ' 1 - 0 , 

RHO  EOUBLS  ()  , .4 

r H e v e e t o F:  m u i s j,  . 4 />  4 7 1 , 4 4 7 l. , A A A 7 1 . A A A 7 1.  * 6 A 6 7 
*************************************************** 
THE  V E C T O R X N O T X S J 

, 0 . 1 2 0 . 1 2 0 . 1 2 0 . 1. 2 0 . 1 2 

THE  V E C TO H X O N E X S ♦ 

0 . 0 5 9 0 1. 9 0.0 5 9 0 1 9 0 . 0 5 9 0 1 9 0.0  5 9 0 1 9 0 - 0 5 9 01" 

**********.*********************  ft  A ft  I-  ft  ft  ft  ft;  * * A A ■'  ft  Ac  A A ft  ft 
THE  CONDITIONAL  MEAN  QUEUE  LENGTH'S  t N TIT  F 

v a e :i  o u s p 1 1 at  s e s a r e * 

0 . 5 3 3 3 3 0 . 5 3 3 3 3 0 ♦ 5 3 3 3 3 0 . 5 3 3 3 3 0 . S 3 3 3 3 


THE  MEAN  QUEUE  LENGTH  IS* 


1 
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A ft  A A A A A * A A A ft  A A A A A A A A ft  A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A 


the: 

vector  IN 

PARTITIONED 

FORM  IS 

G 

xv till  by; 

0 

1.2000E-1 

1 .2000E--1 

1 « 2000E" 

"i 

1 . 2000*2 

“1 

1 .2000*2" 

"1 

t 

5.9(>19E“2 

5.9019*2~2 

5 • 9019*2" 

“p 

5.9019*2 

~2 

5, 9019E 

2 

o 

1 .6437E  2 

1 .6437 *2  "2 

1 • 64  3 7 *2 

.. 

1 . 6437E 

1a 

1 . 6437*2' 

2 

‘7 

O 

3 * 6252*2"“  3 

3 . 6252E" 3 

3 . 6252*2' 

"3 

3 . 625 2 E 

~3 

3 . 6252*2 

"3 

4 

N: 

W 

fT: 

1 

7 . 335 4 E”" 4 

7 * 3354E" 

"4 

7.3354E 

"4 

7.3354*2 

'4 

C’ 

X.J 

1 . 4460e“4 

1 . 4460*2 '“4 

1 . 4460*2" 

4 

1 . 4460E 

“4 

1.4460*2 

'4 

6 

2.7718e"5 

2. 7718E“5 

2 ■ 7718e 

*•  cr 

2*7710*2 

— i:r 

2. 7718*2 

t'j 

7 

4 . 5326*2" 6 

4 . 5826 E~ 6 

4 . 5826*2 

••  / 

U 

4 . 5326*2 

”6 

4.5326*2 

6 

A A A A A A A A A A A A A A A A A A A A A A * A A A * A A A A A :A  A A A A A A A A A A A A A A A A :A  A 


In  the  first  run.  Cl  is  much  larger  than  C.  As  one 
would  expect,  the  conditional  mean  queue  length  is  largest 
immediately  following  long  services  and  gradually  decreases 
until  a minimum  is  reached  at  the  beginning  of  the  long  ser- 
vice. In  the  second  run.  Cl  is  much  smaller  than  C.  Here 
we  note  the  opposite  effect.  The  queue  length  immediately 
following  the  very  short  service  is  smallest,  as  is  to  be 
expected.  In  the  last  case,  C = Cl  and  we  have  an  M/D/1 
queue  with  p = C.  Note  that  in  this  case,  the  mean  queue 


A 
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length  agrees  with  that  given  by  the  classical  formula 

0 + rii^T  ■ 

The  small  scale  implementations  of  the  algorithm, 
selected  for  inclusion  here,  do  not  fully  illustrate  its 
power  to  handle  the  high  orders  of  the  matrices  likely  to 
arise  in  practical  situations.  Examples  with  m as  large  as 
fifteen  were  generated  to  test  our  APL  program  and  runs 
with  m as  high  as  fifty,  using  a well-written  FORTRAN  pro- 
gram, are  entirely  feasible  without  requiring  prohibitively 
expensive  processing  times,  except  in  cases  where  p = rt_B  is 
very  close  to  one.  In  the  latter  cases,  the  practical 
value  of  steady-state  distributions  is  in  fact  questionable. 
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APPENDIX  I 


Proof  of  Theorem  4: 


By  substituting  Equations  (24)  into  (27)  we  obtain 


( 1 ' ) d = 1MV.  kC0(I-B0)_1. 


( < K * ) 


(dd*) 

Since  V— _ . is  a constant,  we  must  first  show  that 

( K k*  ) 

<Cq(I-Bq)  1 is  a left  eigenvector  of  L=L(1).  Multiplying  L 
on  the  left  by  £Cq  ( I-Bq ) ~-*- , we  have 


(2')  kC q ( I - B 0 ) = <C0(I-B0)_1B0  + <C0(I-B0)-1 


(i 


V v-1  ” v_l  -1  1 

' BvG  (I-  l AVG  ) CQ  ‘ 


v*l 


’J 


Rearranging  the  equation  <K  = yields 


_ oo  _ 

(3*)  kC0(I-B0)_1  l BvGV~  = k 

v = l 


Substitution  into  (2')  yields 


r oo 

I 

v=l 


v-1 


(4')  kC0(I-B0)  1L  = kC0(I-B0)  1Bq  + kCq  = kC0(I-Bq)  \ 


Now  since  de=l , it  remains  to  verify  that 


■ ♦ 
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(5')  kk*  = (dd*) <C0(I-B0)e 


(by  multiplying  equation  (1')  on  the  right  by  e) . 

Using  Theorem  3 and  multiplying  the  equation  for  k*  on  the 
left  by  we  have 


(6')  kk*  = 1 + kCq  (I-Bq)  _±e  + k_<Cq  ( I-Bg) 


00 

1 

_v  = l 


Bv  ~ I BVG  + l (v-l)BvG 
v=l  v=2 


1 + I Av 

J V = 1 


- I AvGv_i  + l (V-1)A  G } (I-G+G)  \ , 
v=l  v=2 


and  similarly 


(7’)  dd*  = 1 + d f l Bv  - l BvGV_X  + l (V-1)B 

L V = 1 v=l  v=2 


(I-G+G)  _1m  + d l BvGV_i(j-  l AvGV_i'1  " 
v = 1 v = 1 

(a  + f U - l AVGV-X  + I (v-1)Avg1 

[ [y=l  v = l v=2  J 


(I-G+G)  1p 


The  equation  dL  = d implies 


(8’)  d = d J.  BvGV_1(I-  l AvGV_1)  ^Cod-Bo)'1  , 
v=l  v— 1 


and  therefore  dd*  can  be  written  as 


dd*  = 1 + d [ ByGV  1[l-  l AVGV  ] 1e 
v=l  v=l 


+ ^J1BvgV'1CI-J1AvGV'13“  |C0(I-B0)_; 

[CO  OO  OO  I 00 

l Bv-  l BVGV~  + l (v-l)BvG  + [Av 
v=l  v=l  v=2  J v=l 


- I AyG  + l (v-1) AVG> (I-G+G)  p. 
v=l  v=2  ! 


Comparing  (6')  with  (9'),  we  see  that  (51)  holds  if 


(10')  rKC0(I-B0)_1e]  d l BVGV  X[l-  l AVGV~] J-1  = k. 

V=1  V=1 


By  a straight  forward  substitution,  it  is  easily  verified 

that  d l B Gv  fl-  £ A Gv  is  a left  eigenvector  of  K. 

v=l  v v=l 

Rearranging  the  equations  dL=d  and  £K=^,  yields  equation  (81) 


-rr1 


dl’)  < = kC0(I-B0)  1 l BvGv“1[l-  l AVGV  X] 

v=l  v=l 
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Multiplying  (11')  on  the  right  by  C0(I-BQ)  1 and  comparing 
with  (8')  we  see  that 

(12')  d = 0kCo (I-B0) -1, 

for  some  constant  0.  But  since  de  = 1 , 

( 13 1 ) 0 = QcCq (I-Bq) . 

Upon  substitution  into  (12'),  multiplying  on  the  right  by 

00  oo 

1 BVGV  fl  - l AVGV  1'|  1e  and  using  (11'),  we  finally 

v = l v = l 

obtain 


APPENDIX  II 


Proof  of  Theorem  6 : 


Differentiating  (7)  once  with  respect  to  z yields 


(1")  X ' ( z ) [ z I-A*  ( z ) ] + X(z)  [I-A*’  (z)  ] = x_  7 B,z 

Uk=l  K 


r , k-1 

+ 2*0.  i,kBkz  — 1A0  ‘ 

k=l 


Letting  z tend  to  1- gives 


(2") 


X1  (1)  (I-A)  + X(l)  fl-  l kAkl  = x0  l B,  + l kBk 

k=l  k=l  ”~uk=l 


-1A0' 


We  note  that  I-A  is  singular  but  I-A+ri  is  non- 
singular and  X'(l)n  = (X'(De)j.  Therefore  (2")  becomes 


(3") 


£'(1)  * UoJ^k  + XoJ^Sfc  - xxA0  - X(l)[l-JikA,3 


( I-A+li)  _1  + (X|(l)e)£  . 


Differentiating  twice  in  Equation  (7)  with  respect  to  z 
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t 


84 


gives 


(4")  X" (z) [zI-A* (z) ] + 2X*  (z) [I-A*' (z) ] -X(z)A*"(z) 


= 2xn  l kBkzk_1  + zxn  l k(k-l)B 


,k-l 


k=l 


k=2 


As  z + 1-  we  have 


(5")  X " ( 1 ) [ I-A]  = X ( X ) A* " ( 1 ) - 2X' (1) [I-A* ' (1) ] 


+ 2x0  l kB,  + xn  l k(k-l)Bk, 
k=l  uk=2 


or 


(6") 


X"  (1)  = <X  (1)  [ k(k-l)Ak-2X'  (1)  fl-  [ kAk~) 
k=2  ^ k=l 


+ 2x0  l kBk  + xn  l k (k-1) Bk  ) (I-A+n) 
k=l  k=2 


-1 


+ (X"  (1)  e)  tt  . 


To  determine  X' (l)e,  we  multiply  (7)  on  the  right  by  u(z) 
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(7")  X ' (z)  u(z)  t z — <S  (z)  ] + X (z)  u'  (z)  [z-6 (z) ] 


+ X(z)u(z)  [1-6  ' (z)  ] 


= Xn  l B,  z u(z)  + zxq  [ kB,  zK_1u(z) 
uk=l  k=l  _ 


+ ZXq  l Bj^Z^U1  (Z)  - X^AgU(z)  - ZXj 

k=l 


which  implies 
(8")  X 


^ 03  00 

‘ (Z)^(2)  = z-6  ( z ) i — 0 l BkzkH(z)  + z*0  ^ 

k=l  k=] 


+ ZXg  l Bj^Z^u'  (z)  - J^AqUCz)  - ZXj 

k=l 


- X ( z)  u ( z)  [ 1-5  ’ ( z)  ] > - X ( z)  u ' ( z ) . 


Letting  z -+  1 yields 


(9")  X' ( 1) e = i 


0 l B ® + *0  l kBk£  + *o  £ 1 


'k=l 


k=l 


k=l 


- x^Age  - x^Ag  u ' (1) 


- X(l)e(l-p)] 


AgU'  (z)  , 

kB^z^'^u (z) 
AgU' (z) 

kH.'  (1) 

- X (l)u'  (1) 


Note  that  the  term  in  braces  equals 
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OO  00  00 

<10">  *0  l Bk®  + *o  l kBk—  + *o  l -Bk  C(I_A+n) 

k-1  k — 1 k=l 

-XiAQe  - x-jAq  [(I-A+n)  _1£-pej  - X(l)e(l-p) 

OO 

= (1-p)  fx0  (e-B0e)  -Xj^AQeJ  + xQ  £ kBke 

k=l 

00 

+ fx0  l Bk-X1A0] (I-A+n)_16  - X (1) e (1-p) 
k=l 

oo 

= xQ  l kBke  + [x  (1)  - (l-xQe)  tT)b  - X(l)e  + X(l)ep 
k=l  J 


x0  l kBke  + X(l)6.  - X (1)  e = 0, 
k=l 


by  multiplying  (3")  on  the  right  by  e and  simplifying. 
Therefore,  we  may  use  L'Hopital's  rule  on  Formula  (8")  to 
get 

(11")  X'  (z)u(z)  = 1-J61(Z)  dXg  l kBkzk-1u(z) 

l k=1 

oo  00 

+ 2xn  l Bkzku' (z)  + zxn  y k (k-1) Bkzk_2u ( z) 

~Uk=l  “ “°k=2  k 


W OO  - 

+ 2zx0  l kBkzk-1u' (z)  + zxa  l Bkz  u"(z) 
k=l  k=l  “ 


i 


- 2x1AqU'  (z)  - zx^j^AqU"  (z)  + X (z)u  (z)  6"  (z) 


- X'(z)u(z)  - 2X(z)u'(z), 


which  finally  yields 
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(12") 


X’ (1) e = 


® = | 2— 0.  kBk—  + 2— 0 ^ Bk— ' 

I k=l  k=l 


(1) 


OO  00 

+ Xq  [ k(k-l)Bke  + 2Xq  l kBku’ (1) 
k= 2 k-! 


+ *0  ^ BkH" t1) 
k=l 

+ X (1) e<5" (1)1 


ZXj^AqU' (1)  - ^AqU"  (1) 

X(l)u'  (1)  . 


In  order  to  evaluate  X"(l)e,  we  differentiate  (7") 
with  respect  to  z to  get 


(13")  X"  (z)u(z)  [z-6  (z)  ] + 2X ' (z)u'(z)[z-6(z)] 


+ 2X'  (z)u(z)  [1-6*  (z) ] + X(z)u" (z)  [z-6 (z)  ] 


+ 2X(z)u'  (z)  [1— 6(z)  ] -X(z)u(z)6"(z) 


00  OO 

= 2Xq  l kBkzk  1u(z)  + 2xq  l Bkzku' (z) 


OO  00 

+ zx0  l k(k-l)B  zk-^u(z)  + 2zx0  £ kBkz^-^u' 
k=2  K k=l 


+ ZXq  l BkZku"(z)  - 2x1AqU' (Z)  - ZX^AqU" ( z ) 


which,  upon  setting  z=l,  simplifies  to 
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(14")  X " ( 1 ) e = 


£ <x(l)e6"(l)  - 2X ' (l)e(l-p) 


- 4>X(l)u' (1)  (1-p)  + 2x0  l kB,e 

k=l 


+ 2x0  l Bku' (1)  + xQ  l k (k-1) B,  e 
k=l  _Uk=2 


+ 2xq  l kBku'(l)  + x0  l Bku"(l)  - 2x1A0u'(l) 
k=l  k=l 


- Xj^AqU"  (1) 


- 2X ' (l)u'  (1)  - X(l)u"  (1)  . 


In  order  to  show  that  the  term  in  braces  equals  zero,  note 
that  from  (10")  , 


(15")  xQ  l kBke  + xQ  l Bku'  (1)  = XjA0e  + ^AqU  ' (1) 
k=l  k=l 


+ X (1)  e (1-p)  -xQ  l Bke 


k=l 


Upon  substitution  of  (10")  and  (15")  into  the  braces  and 

oo 

noting  that  XjA0e  - xQ  l Bke  = - Xq (e-BQe)  = 0 , we 

k-1 

see  that  the  term  in  braces  reduces  to  0.  Therefore,  we  can 
again  apply  L'Hopital's  rule  to  get 
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(16")  X"  (z) u (z)  = 


3(1-6'  ( z) ) 1 - 


3X ' (z)u(z) 6" (z) 


+ 3X(z)u'(z),s"(z)  + X(z)u(z)  a'"  (z) 


00  , 00  , _ 1 

+ 3xq  l k (k-1)  Bjczku  (z)  + 6Xq  £ kB^z*  u*  (z) 


k=2 


k=l 


00  00 

+ 3xq  l Bkzku"(z)  + zXq  l k(k-l)  (k-2)Bkzk-3u(z) 


k=l 


k=3 


00  w 

+ 3zxq  l k (k-1)  Bkzk_2u ' (z)  + 3zxq  £ kBkz^_1u"(z) 
k=2  k=l 


+ ZXq  l B,zku’"  (z)  - 3xi  Anu"  (z) 
k=l 


- ZX-^AqU"'  (z) 


\ - | X' (z)u' (z)  - X ( z ) u " (z) . 


Letting  z -*■  1-  yields  Equation  (44)  and  completes  the  proof. 

The  higher  factorial  moments  of  the  queue  length 
may  in  principle  be  computed  in  the  same  manner  but  the 
formulas  become  uninspiringly  complicated  and  will  not  be 
shown  here. 
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